Paul Terzian

Eric Olo Ndela

Clovis Galiez

Julien Lossouarn

Rubén Enrique

Pérez Bucio

Robin Mom

Ariane Toussaint

Marie-Agnès Petit
email: marie-agnes.petit@inrae.fr

François Enault
email: francois.enault@uca.fr

Rub

Enrique P Érez Bucio

Marie-Agn Ès Petit

PHROG: families of prokaryotic virus proteins clustered using remote homology

come

INTRODUCTION

Viruses are key players in most ecosystems as they actively participate to the regulation of microbial communities [START_REF] Breitbart | Here a virus, there a virus, everywhere the same virus?[END_REF] and are important vectors for horizontal gene transfer [START_REF] Suttle | Marine viruses --major players in the global ecosystem[END_REF]. Viruses infecting prokaryotes out-number eukaryotic viruses in some ecosystems [START_REF] Reyes | Viruses in the faecal microbiota of monozygotic twins and their mothers[END_REF] and represent the great majority of the viruses found in viral metagenomes [START_REF] Gregory | Marine DNA viral macro-and microdiversity from pole to pole[END_REF]. Viral genomes and metagenomes highlight (i) the recurrent presence of virus groups, especially in similar ecosystems, recently termed 'virus operational taxonomic units' (vO-TUs; 5), (ii) the large number of different vOTUs found in most ecosystems [START_REF] Gregory | Marine DNA viral macro-and microdiversity from pole to pole[END_REF] and (iii) the dominance in most vO-TUs of genes not similar to any known genes. The growing amount of viral sequences produced nowadays, especially from metagenomes, calls for a need in developing resources to help assigning functions to viral proteins in order to improve functional and comparative analyses. As determining information for newly identified genes in sequences relies on finding an annotated homologous sequence using similarity searches, a set of well-annotated proteins that can be used for future work is particularly important. In order to build such a reference set of proteins, a classical way is to organize these proteins into homologous groups and to annotate these groups. Several methods have been used to cluster viral proteins into either homologous or orthologous groups: (i) an approach based on the identification of genome-specific best hits that are joined to form clusters of orthologs has been used for the pVOGs database (Prokaryotic Virus Orthologous Groups; [START_REF] Grazziotin | Prokaryotic Virus Orthologous Groups (pVOGs): a resource for comparative genomics and protein family annotation[END_REF]), (ii) a similar approach based on best-hit triangles has been implemented recently in eggNOG (evolutionary genealogy of genes: non-supervised orthologous groups; [START_REF] Huerta-Cepas | eggNOG 5.0: a hierarchical, functionally and phylogenetically annotated orthology resource based on 5090 organisms and 2502 viruses[END_REF]), which integrates an additional step of in-paralogs detection and the identification of fused genes and (iii) clustering-based approaches have also been used to compute groups of homologous viral proteins [START_REF] Li | OrthoMCL: identification of ortholog groups for eukaryotic genomes[END_REF][START_REF] Roux | Viral dark matter and virus-host interactions resolved from publicly available microbial genomes[END_REF][START_REF] Soares | GIPSy: Genomic island prediction software[END_REF]. Despite using different clustering strategies, all these methods rely on similarity search results generated by BLAST [START_REF] Altschul | Basic local alignment search tool[END_REF] or faster equivalent such as MMseqs [START_REF] Steinegger | MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets[END_REF] or DIAMOND [START_REF] Buchfink | Fast and sensitive protein alignment using DIAMOND[END_REF], and none integrate remote homology detection. Yet, viruses are known to have distant evolutionary relationships, resulting in distant sequence similarities not always captured by sequence comparison tools. This can be illustrated by the Microviridae family, whose members all encode homologous major capsid and replication initiation proteins in their ∼5 kb genomes. As the origin of the family is ancient, homologs are difficult to detect using sequence similarity search tool such as BLAST. For instance, for two distantly related Microviridae, Spiroplasma phage 4 and Enterobacteria phage phiX174 belonging respectively to the Gokushovirinae and Bullavirinae sub-families respectively, the best BLASTp hit between their capsid proteins only exhibit a bit-score of 38.5 (corresponding to an E-value of 0.14 on a ∼1 million protein database) and an alignment coverage below 30% for the two proteins. These values are much below any reasonable thresholds on the bit-score, E-value or coverage to detect homology. Thus, no clustering methods will group these two proteins together even though they are distant yet true homologs. Recent developments in remote homology detection software [START_REF] Remmert | HHblits: lightning-fast iterative protein sequence searching by HMM-HMM alignment[END_REF][START_REF] Steinegger | HH-suite3 for fast remote homology detection and deep protein annotation[END_REF] coupled to the great wealth of genomic sequences could help tackle the problem of distant homology of viral proteins, yet to our knowledge, no clustering method uses remote homology detection in an automatic way for viral proteins.

Here, we propose such a procedure and clustered viral proteins into homologous groups, in two steps: (i) proteins are first gathered based on similarity search results (score >30 and coverage >80%) and (ii) HMM profiles generated for each protein cluster are then compared to each other and grouped (coverage >60%, probability >90 %) into superclusters hereafter termed PHROGs (Prokaryotic Virus Remote Homologous Groups). These PHROGs were then annotated using annotation transfer, coupled with a careful manual curation. An interface has been developed, allowing to browse the data either by PHROG or by viral genome.

MATERIALS AND METHODS

Two datasets of archaeal and bacterial viruses

First, all known sequences of viruses infecting Bacteria or Archaea were assembled in a dataset. To this end, 2318 reference sequences of viruses infecting prokaryotes were retrieved in RefSeqVirus genomes (as of April 2018). Viruses included in the pVOGs database [START_REF] Grazziotin | Prokaryotic Virus Orthologous Groups (pVOGs): a resource for comparative genomics and protein family annotation[END_REF] and complete virus genomes from GenBank were also downloaded. As one virus can be present in these three databases with different identifiers (RefSeq, pVOGs and GenBank), all these viruses were compared to each other using BLASTn and viruses almost identical were removed (>99.9 identity percent and coverage >97%). The resulting 4975 completely sequenced viruses (2315 from RefSeqVirus, 686 from pVOGs and 1986 from GenBank) will be refered as RefVirus from here on. To this dataset, 12 498 previously published curated viral sequences derived from cultivated microbial isolates were added [START_REF] Roux | Viral dark matter and virus-host interactions resolved from publicly available microbial genomes[END_REF]. This last dataset is composed of both integrated proviruses (>10 kb) and circular episomes, and is here termed ProVirus. These (pro)viruses were found in 5492 microbial genomes. In total, 496 859 proteins from complete viral genomes (i.e. RefVirus) and 442 005 proteins from (pro)viruses (i.e. ProVirus) were collected.

Generating the viral homologous groups or PHROGs

Protein clustering using similarity searches (Figure 1, Panel B1 and B2). The 938 864 proteins were compared to each other using MMseqs [START_REF] Steinegger | MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets[END_REF]. To be further considered, a protein pair should have (i) at least a local alignment with a bit-score >30 and (ii) >80% of the residues of each protein should be involved in at least one alignment found between the two proteins (i.e. coverage >80%). Using each protein pair and the lowest E-value found in a local alignment between two proteins, the proteins were clustered with MCL (inflation 2.0; ( 16)).

Grouping the protein clusters using remote homology detection (Figure 1, Panel B3 and B4). For each of the 63 673 clusters containing at least two proteins, a multiple alignment was built using ClustalOmega (17; 5 guide-tree/HMM iterations) and an HMM profile was computed for each alignment using hhmake of the HHsuite toolkit (version 2.0.16; (18)). All these profiles were then compared to each other using the hhsearch command. The 85 024 singletons were also compared to the 63 673 cluster profiles. To be further considered, a cluster pair should have a hit (i) with a probability >90%, (ii) that involves at least 60% of the two HMM profiles (same thresholds when comparing singletons to clusters). Based on these rules, singletons and clusters were then clustered with MCL (inflation 2.0; [START_REF] Enright | An efficient algorithm for large-scale detection of protein families[END_REF]). This resulted in placing 868 340 of the initial 938 864 protein sequences into 38 880 'super-clusters' containing at least two proteins, hereafter named PHROGs. Only 7.5% of the protein dataset remained as singleton, or ORFan (70 524 proteins). Multiple sequence alignments and HMM profiles were computed for all the 38 880 PHROGs containing at least two proteins using ClustalOmega [START_REF] Sievers | Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega[END_REF] and HHsuite [START_REF]Protein homology detection by HMM-HMM comparison[END_REF]. Each multiple alignment was processed in order to consider as gaps the columns made of >50% of gaps (saved in a2m format). HMM profiles were generated for each of the masked alignment and these profiles were compared to each other using hhsearch [START_REF]Protein homology detection by HMM-HMM comparison[END_REF]. Results are accessible through the website and all files are downloadable as zipped archives.

Comparing PHROGs to the standard clustering method

To estimate the performance of our procedure, the same set of proteins were clustered using a classical approach. Sequence similarities were searched for all protein pairs using MMseqs and clustered using Markov clustering algorithm (here the MCL software) [START_REF] Steinegger | MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets[END_REF][START_REF] Enright | An efficient algorithm for large-scale detection of protein families[END_REF]. This clustering procedure is similar to the first step of the PHROG clustering described above, except that a less stringent coverage threshold was used (coverage ≥50 % for the two proteins). Two complementary metrics were used to analyze the differences between the two methods: the number of proteins clustered and the size of the clusters (Figure 3), and the identity percent and coverage between proteins grouped in each cluster (Figure 4). To this end, multiple alignments of clusters generated with the classical approach were also computed using ClustalOmega. For all protein pairs inside these multiple alignments and the PHROGs, the percent identity and the coverage were calculated as being, respectively, (i) the number of identical amino acids in the two aligned proteins divided by the length of the smallest of the two proteins and (ii) the proportion of amino acids of one protein that is aligned to any amino acid (not to a gap) of the other protein. Subsamples of 1000 values were taken to draw each boxplot of Figure 4.

Coverage threshold values chosen for the two steps of the clustering procedure are critical parameters, since relaxed values come with the risk of aggregating heteroge-neous functions (due to the presence of protein fusions), and stringent values limit cluster sizes. Based on a subset of 50 000 randomly extracted proteins, three thresholds on coverage (80%, 50%, no coverage threshold, named respectively '1-step C80', '1-step C50','1-step C0') were tested for the first step of the clustering procedure (the 50 threshold corresponding to the standard procedure described above). Then, results from the 80 % threshold on coverage were further considered and used as a basis for the second step of the clustering and 5 thresholds were tested for the coverage between HMM pairs (80%, 60%, 40%, 20% and no threshold, named respectively '2-step C80 C80', '2-step C80 C60', '2-step C80 C40', '2-step C80 C20', '2-step C80 C0'), the combination of C80 and C60 being the ones retained for building PHROGs.

Functional annotation of PHROGs

NCBI protein annotations (RefVirus) were first automatically curated (no upper case except in gene names, correction of typos, etc.). Phages for which all protein annotations was only an uninformative list ('hypothetical protein GP1', 'hypothetical protein GP2', etc.) were considered as unannotated. These protein annotations were first combined to determine an annotation for PHROGs. To refine this annotation, the 38 880 PHROGs were compared to different databases. PHROG profiles were compared to Pfam domains (version of jan 2018; 19) and UNICLUST [START_REF] Mirdita | Uniclust databases of clustered and deeply annotated protein sequences and alignments[END_REF] and individual viral protein were compared to proteins in KEGG Orthologous groups (KOs; version of jan 2018) using MMseq (bit-score>50, coverage >50%). Manual curation of the collected annotations and similarities allowed to extract a single annotation per PHROG.

Colocalization of genes according to their functions

For each pair of functional annotations, the number of their respective proteins that are neighbors on the genomes were determined. Neighbors were here defined as the two proteins encoded right before and after the protein considered, and only if they are on the same strand. As some annotations gather many more proteins than others, a colocalization score between two annotations was computed. The hypergeometric formula was first used to calculate the probability P that two annotations (noted here annotation A and annotation B) that each have n A and n B proteins (among the whole dataset of n proteins), would have k pairs of proteins that are neighbors by chance.

P(X ≥ k) = mi n(n a ,n b ) i =k C i n a C n b -i n-n a C n b n
This probability represents the statistical significance of an observed number of neighbor protein pairs for two annotations. This probability was also computed considering each PHROG and not each annotation and PHROG neighbors are reported on the Web site. To draw a graph of the most significant relationships, all P values were then corrected by the total number of annotations pairs T and finally converted into a significance score : S(A,B) =log(P×T) (as in 21). Annotation pairs with significance scores >1 are considered as significantly colocalized (scores are capped to 1000). A network of neighboring annotations was generated using the igraph library in R, in which annotations were joined by an edge with a weight equal to their significance score. For visualizing purposes, only the 112 more frequently retrieved annotations (used for >650 proteins) were considered. When neighbor genes have the same annotations, these annotations are drawn as squares. Annotations were placed using the Fruchterman-Reingold layout algorithm. Annotations were gathered into modules using community structure detection based on edge betweenness.

RESULTS

Description of the dataset

All prokaryote-infecting viral sequences (complete genomes or not) from existing sequence databases were compiled. Overall, to the 4975 viruses infecting a prokaryote from the RefSeq, pVOGs and GenBank databases were added the 12 498 high-confidence viral sequences identified by the VirSorter tool [START_REF] Roux | Viral dark matter and virus-host interactions resolved from publicly available microbial genomes[END_REF][START_REF] Roux | VirSorter: mining viral signal from microbial genomic data[END_REF] as proviruses or circular viral genomes in microbial genomes. The viruses considered here were mostly dsDNA viruses and infected 35 different bacterial and archaeal classes and 410 genera. Yet, viruses infecting bacteria belonging to the Proteobacteria class represented as much as 63% of the 17 473 total viral genomes or contigs, the Escherichia genus alone representing as much as 29% of this set (Figure 2).

Analysis of the clustering characteristics and comparison to a standard clustering procedure

Viral proteins were clustered using either our two-step procedure or a more classical approach (see Materials and Methods). The main improvement provided by the PHROGs strategy was a reduction by nearly two-fold in the number of clusters (at least 2 proteins) compared to the standard procedure (38 880 and 57 073, respectively) and a concomitant increase in the average size of the protein clusters (22.3 and 15.2, respectively). Unsurprisingly, the number of proteins left as singletons was comparable for the PHROG and standard procedures (70 524 and 67 430, respectively) and represented <8% of the 938 864 proteins. The largest clusters were composed of 5879 and 2267 proteins, respectively, with 19 PHROGs being >2000, whereas only 2 clusters made with the classical approach had such a large size. All along, highly populated protein families were more frequent among PHROGs than traditional clusters (see Figure 3).

Considering the classical clustering approach (Figure 4), the average identity percent inside clusters tended to decrease with the cluster size, from 79.7% for small clusters (<7 proteins) to 56.6% for large ones (>2000 proteins). This trend was dramatically accentuated for the PHROG procedure, as the average identity went down from 79.8 to 27.2%. The average coverage remained high for all cluster sizes for the classical approach (98.3% in average) whereas it dropped for the PHROG procedure from 99% to 80.7% in average, for small and large clusters, respectively (Figure 4). This shows that more distant homologs are gath- ered with the PHROG procedure. We verified that these large clusters maintained a good consistency, with an acceptable average coverage: over 95.8% of the proteins inside a large PHROG covered >50% of any other protein of the PHROG. As expected, the percent identity were quite low within these large clusters, with 66.8% of the protein pairs having an amino acid percent identity <30%. Coming back to the Microviridae protein example mentioned in the Introduction, capsids of viruses of the two known sub-families, Gokushovirinae and Bullavirinae, were grouped into two separate clusters by the first step of our method, as well as by the classical method. Indeed, using pairwise sequence alignment (MMseqs), all Gokushovirinae capsid protein sequences were similar to each other but no similarities were detected between these and Bullavirinae capsid sequences. An HMM profile was built for these two groups of capsids and the two HMMs were clearly identified as similar by HHsearch (probability of 99.5, E-value of 8.4 × 10 -18 , score of 176.8, covering >400 positions of the two HMMs). Thus, the two clusters of capsid proteins were gathered by the second step of our method, resulting in PHROG 514 that contains capsid protein for all members of this phage family.

We further studied the effect of coverage thresholds on clusters and found that all procedures based on only one clustering step gave comparable results, regardless of the coverage threshold (Supplementary Figure S1A). This parameter however had a major impact on the second clustering step (Supplementary Figure S1A). Indeed, for '1-step C80', '1-step C50' and '1-step C0', the largest clusters were of 71, 77 and 90 proteins respectively, whereas the 2-step procedures built clusters composed of 151, 164, 269, 465 and 343 proteins when coverage thresholds of the second step were decreased from 80% to 0, by 20% steps, using a coverage of 80% in the first step. Considering the average identity percent for protein pairs inside clusters, the use of HMM comparisons allowed to group proteins with significantly more remote homology (Supplementary Figure S1B). Yet, procedures involving thresholds lower than 50 % on the first or second step ('1-step C0', '2-step C80 C40', '2-step C80 C20' and '2-step C80 C0') of course resulted in the non negligible presence of protein pairs that cover each other <50% inside clusters (Supplementary Figure S1C). Thus, the combination of a threshold of 80% and 60% on coverage for the two steps of the clustering allows to detect remote homologies (<30 id%) and to build large clusters without grouping proteins that do not cover themselves significantly.

The whole protein dataset, when compared by a sequence-to-sequence comparison tool (MMseqs, bit-score>50), resulted in 238.3 millions of protein pairs. Considering only the 195.7 million pairs with a coverage >50%, the standard method built 57 073 clusters (178.5 millions of protein pairs inside these clusters). The 38 880 PHROGs created by our clustering strategy resulted in 463.1 millions of protein pairs inside these PHROGs, so over twice as much the number of pairs identified by BLAST-like searches. Using the multiple alignments of each cluster, (i) the identity percent between two proteins is the number of amino acids that are identical in the two aligned proteins divided by the length of the smallest of the two proteins, and (ii) the coverage is the proportion of amino acids of one protein that is aligned to any amino acid (not to a gap) of the other protein. Subsamples of 1000 values where taken to draw each boxplot.

Functional annotation of PHROGs

To determine a putative function for each PHROG, protein annotation from NCBI complete viral genomes (Re-fVirus) inside each PHROG were confronted, harmonized and transferred. At the onset, as many as 20 872 different annotations were found for the 496 859 proteins from Re-fVirus. For example, PHROG 2 is made of proteins annotated under 107 different designations such as 'terminase large subunit', 'terminase', 'terminase DNA packaging enzyme large subunit', 'phage terminase large subunit', 'large terminase protein', 'large terminase subunit', 'TerL', 'terminase, large subunit', 'large subunit terminase', 'large terminase', (etc. . . ) as well as some 'hypothetical protein'. This PHROG was then annotated as 'terminase large subunit'. To further improve these annotations, different sources of information were collected and manually inspected to refine each PHROG annotation: (i) 11 810 PHROGs contain at least one Pfam domain, (ii) 6710 PHROGs were attributed at least one GO term via a similarity to UNI-CLUST proteins, (iii) 1788 PHROGs were similar to at least one KEGG Ortholog group (KO). Presence of proteins annotated unambiguously in the ACLAME database [START_REF] Leplae | ACLAME: A CLAssification of Mobile genetic Elements, update 2010[END_REF] were also retrieved. Finally, 5108 of the 38 880 PHROGs were functionally annotated using a set of 705 annotations. Nine functional categories were then defined to cover these 705 annotations : 'head and packaging', 'connector', 'tail', 'DNA, RNA and nucleotide metabolism', 'integra-tion and excision', 'lysis', 'transcription regulation', 'moron, auxiliary metabolic gene and host takeover', 'other' (Table 1 and Figure 5). Even though the number of annotated PHROG seems low (5108 of the 38 880 PHROGs), large PHROGs are more often annotated than small ones and thus, these 5108 annotated PHROGs contain 475 493 proteins (50.6% of the total protein dataset, including singletons) (Table 1). These annotated PHROGs contain as many as 68 478 'hypothetical protein' from RefVirus that can now be annotated with their PHROG's annotation. On the contrary, 11 660 RefVirus proteins initially considered as annotated have lost there annotation (i.e. are in unannotated PHROGs; Figure 1C). This is due to PHROGs (i) containing only uninformative protein annotations (e.g. protein annotations such as 'prophage Lp1 protein 30', 'DUF1642 domain containing protein', 'similar to P2 orf80', etc. . . ) or (ii) containing protein with non congruent annotations and in both cases, information collected using KEGG, UNI-CLUST or Pfam did not help to define an annotation. For example, PHROG 358 contains RefVirus proteins annotated as 'virion structural protein', 'virulence associated protein', 'baseplate protein', 'neck protein', 'tail protein' (etc.), has no similarity with KEGG or UNICLUST proteins, is only similar to a Pfam domain of unknown function (DUF4815) and is similar to five unannotated PHROGs and to one 'virion structural protein' PHROG but on a very small portion. Thus, this PHROG was attributed to 'unknown function'.

Further gaining some functional clues using (often short) similarities between PHROGs

We then performed iterative HMM alignments, with relaxed parameters, to extend annotation possibilities. For this, columns of the multiple sequence alignments with >50% gaps were considered as insertion states, allowing to get rid of spurious ends and insertions present in only a handful of proteins. A total of 26 177 new connections were detected between 12 791 PHROGs (E-value < 0.001 and no threshold on the coverage), and in 25% of cases, non annotated PHROGs were joined to annotated ones. These PHROGs were not clustered together in the first place, for several reasons. First, similarities had to cover at least 60% of the two initial protein clusters (Figure 1B(iii)) to build PHROGs and thus, local similarities existing between protein clusters were considered insufficient for aggregation. Moreover, MCL clustering (Figure 1B(iii)) implies the presence of more edges within members of PHROGs than between different PHROGs, but links between clusters (and thus PHROGs) exist. Finally, PHROGs are larger and more informative than the protein clusters they are built upon and thus homologies even more distant can be detected when comparing PHROGs than when comparing initial clusters. In any case, similarities between PHROGs can provide hints, especially for unannotated PHROGs similar to annotated ones. Among annotated PHROGs, 3409 different PHROGs were similar to at least one other annotated PHROG and 52.5% of the 10 752 similarities involved PHROGs with the same annotation (83.4% with the same functional category). Moreover, among the 824 similarities for which the two PHROGs reciprocally covered each other by >60% of their length, 64.7% were involving PHROGs with the same annotation. This suggests that information provided by these more distant links could be of relevance for exploring the function of unannotated PHROGs. Indeed, among PHROGs of 'unknown function', for which gaining some information is essential, 3037 (48 704 proteins) were linked to annotated PHROGs (8179 links). These information are indicated in the lower right panel of each PHROG individual page on the Web page.

Further gaining some functional clues using colocalization of annotations

To further help in the annotation of PHROGs, a score of colocalization was computed and was significant for 1731 annotation pairs. For example, the 5122 and 13 119 genes in PHROGs annotated as 'head scaffolding protein' and 'major head protein' were adjacent 4690 times on the different viral genomes and the link had the maximum weight of 1000. Colocalized annotation terms often involved annotations from the same functional category (Figure 6). Some of them corresponded to known functional interactions, such as 'integrase' and 'excisionase' [START_REF] Cho | Interactions between integrase and excisionase in the phage lambda excisive nucleoprotein complex[END_REF], 'SbcC-like subunit of palindrome specific endonuclease' and 'SbcD-like subunit of palindrome specific endonuclease' (25), 'RecTlike ssDNA annealing protein' and 'exonuclease' (for a review see 26) or 'Sak4-like ssDNA annealing protein' and 'single strand DNA binding protein' [START_REF] Hutinet | Sak4 of Phage HK620 Is a RecA remote homolog with single-strand annealing activity stimulated by its cognate SSB protein[END_REF]. Other, less expected colocalizations, such as 'UvsX-like recombinase' and '2OG-Fe(II) oxygenase' or 'gam-like host nuclease inhibitor' and 'Kil protein for bacterial septation inhibition', might suggest new avenues of investigations.

DISCUSSION

No ORFans anymore?

The PHROG database includes 17 473 genomes (or genome fragments) of viruses infecting a broad range of prokaryotic hosts (410 prokaryotic genera). The great majority of the 938 864 proteins encoded within these genomes could be assigned to a PHROG, and only 7.5% of them remained as singletons. In accordance with previous studies, these 70 524 ORFans were shorter than the other proteins, with an average length of 125.3 amino acids compared to 216.1 for all Figure 5. For each annotation term, the number of PHROGs with this annotation and the number of proteins in these PHROGs, each term being colored according to its functional category.

proteins. It should be noted that our cutoffs for PHROG inclusion were stringent, in particular in terms of alignment coverage (>60%), and lowering this threshold would have allowed to affiliate 23 937 more of these 70 524 singletons. These singletons similar to a PHROG might correspond to protein fusions or to proteins for which the gene prediction software wrongly predicted the beginning or end. For the remaining 46 587 proteins (4.96% of the whole dataset), a non-negligible fraction of these ORFans might not be bona fide genes but false predictions, sometimes due to errors in genomes, as single nucleotide insertions or deletions lead to spurious gene calling. Moreover, for the viruses found in RefSeq and GenBank (RefVirus), various gene prediction programs were used to annotate their genomes over the years and these programs display small differences in the genes they predict for a given genome. Thus, <5% of the proteins can be considered as real ORFans, a value much lower than what is usually described. Their great number in genomes are often described as being a specificity of viruses, with approximately 30% of a phage proteins being described to be such ORFans [START_REF] Frost | Mobile genetic elements: the agents of open source evolution[END_REF]. As already mentioned, the viruses considered here were not randomly sampled, as for example viruses infecting Escherichia represent 29% of the current dataset. Yet, 25 viral families infecting 410 prokarotic genera were present and we can thus assume that ORFans are rare in viruses and that their over-representation in newly described viruses is only due to a lack of closely related reference in databases. This proportion of ORFans much smaller than previously thought helps to draw this picture on the genetic diversity in prokaryotic viruses: viral protein families are numerous and internally genetically diverse, yet these families are conserved in related viruses.

The two steps of the clustering procedure and the search for remote homology

The standard procedure usually implemented to build protein clusters consists in comparing all proteins using a sequence similarity search tool and then applying a Markov clustering algorithm to define groups [START_REF] Roux | Viral dark matter and virus-host interactions resolved from publicly available microbial genomes[END_REF][START_REF] Enright | An efficient algorithm for large-scale detection of protein families[END_REF]. As many as Annotations for which genes were colocalized with genes having the same annotation are drawn as squares. Among these 112, the 21 annotations not significantly colocalized to any other are displayed on the right. 47% more clusters were defined by such a classically used procedure, compared to PHROG, which leads to an overestimation of the number of viral protein families and thus of viral genetic diversity. The clusters delineated by the standard procedure being of a smaller size than PHROGs, the average protein coverage and identity percent inside clusters were greater for this classical procedure. Most protein pairs inside large PHROGs (>2000 proteins) had an identity percent below 30%, a threshold described as the twilight zone under which protein pairs are historically considered to have non-similar structure [START_REF] Rost | Twilight zone of protein sequence alignments[END_REF]. Yet, the coverage percentages inside PHROGs were high, even for large PHROGs, and combined to the threshold of 90% on the probability of HHsearch, indicate that protein inside these PHROGs are bona fide homologs. The low amino acid percentages likely reflect the high diversity of each viral protein family and not the fact that the grouped protein have a non similar structure. Indeed, studies on particular viral protein families have shown that some proteins sharing as little as 10% amino acid identity still belonged to the same functional category [START_REF] Lopes | Automated classification of tailed bacteriophages according to their neck organization[END_REF]. Thus, PHROGs are larger than traditional clusters and the coverage between proteins inside each PHROG is still good. To explain this result, we believe that the first step of the PHROG procedure, although similar to the standard clustering procedure, builds very ho-mogeneous protein clusters as the very stringent threshold (80% of coverage for the two proteins) prevents transitivity issues. Then, the second step based on the detection of distant homologies, allows to securely connect remote homologs.

More viral proteins are now annotated thanks to the annotation of PHROGs

About 13.8% of the RefVirus proteins that were 'hypothetical protein' are now in annotated PHROGs. Even more proteins are newly annotated thanks to PHROGs, as a not negligible (but difficult to estimate) fraction of proteins considered as annotated in RefVirus are not informative (e.g. 'similar to P2 orf80'). These proteins were not detected by our script as 'hypothetical protein' even though no real information on their putative function is known. Another positive improvement (also difficult to quantify) is the correction of existing annotations, that could be misleading, partial or wrong. For example, 3 proteins in PHROG 2, annotated as 'terminase large subunit', were initially annotated as 'type I secretion target domain-containing protein' or 'Mu portal protein gp28 TerL' and annotating these 3 proteins as 'terminase large subunit' feels more informative and legitimate. Furthermore, 263 777 proteins of the ProVirus set are now annotated in a uniform way.

Hints for unannotated PHROGs similar to annotated ones

At present, PHROG annotations are not transferred to unannotated PHROGs that are collected when an additional HHsearch comparison of PHROG profiles is run (no threshold on coverage). However, these results are shown on each PHROG web page and should provide valuable hints. Indeed, more than half of the time, annotated PHROGs similar to other annotated PHROGs have the same annotations. The two most frequent cases are PHROGs 'tail length tape measure protein' similar to PHROGs with the same annotation (1044 different PHROG pairs) and 'tail fiber protein' PHROGs similar to other 'tail fiber protein' PHROGs (929 cases). These two cases are coherent as (i) 'tape measure' have different lengths by nature and (ii) 'tail fiber proteins' being the point of contact with bacteria, they are prone to natural selection, and therefore highly variable in essence. For PHROGs similar to a PHROG with a different annotation, the annotations are often coherent, one being more specific than the other such as 'tail protein' PHROGs similar to 'tail fiber protein' PHROGs (433 cases). Some PHROGs from different categories are also linked because they share a functional domain, such as 'tail length tape measure protein' ('tail' category) and 'endolysin' ('lysis' category) (65 cases). These similarity results should help predict in the future the function of the 3037 unannotated PHROGs (48 704 proteins) that are similar to an annotated PHROG.

CONCLUSION

Due to the ancient origin of viruses, their gene families have a long evolutionary history and often encompass distant homologs not detected by standard sequence comparison tools such as BLAST. To be able to identify these distant homologs, our clustering strategy involved the use of HMM comparison tools as these are much more sensitive. The clusters built here proved to be larger while remaining cohesive, leading to an increase in annotated proteins. In addition, a special effort was made to standardize these annotations. Annotations of PHROG families were performed based on Refseq and several other databases, and manually curated by experts, adding a real value to the PHROG database. These annotations will be updated in the future, as new phage functions are discovered, or when experts provide 'Suggestions' through the dedicated page of the web site. Hopefully, PHROG will constitute a useful tool for scientists to better characterize new viral sequences, especially the millions of proteins derived from metagenomes, bringing useful insights on the nature of viruses infecting microbes.

Figure 1 .

 1 Figure 1.Overview figure of the dataset, the clustering procedure, results and website main page. (A) Protein sets of reference viruses from the NCBI (RefVirus) and (pro)viruses detected using VirSorter (ProVirus) were collected. It should be noted that Provirus proteins (green) were not initially annotated. (B) The four steps of the clustering procedure: (i) the protein network built from pairwise sequence similarities, each dot/vertex representing a protein (green for ProVirus proteins, red for annotated RefVirus proteins and pink for unannotated RefVirus), linked by edges if the two proteins are similar. (ii) Protein clusters are identified by applying MCL on to this network, and clusters are depicted as gray circle. (iii) Clusters (and singletons) are compared using protein profiles and edges are drawn for pairs of protein clusters that are similar. (iv) This network of clusters is here again clustered into PHROGs, depicted as dark brown and blue for unannotated and annotated PHROGs. For example, PHROG2 is made of two protein clusters and one singleton. (C) Description of the number of annotated and unannotated PHROGs and singletons, with the number and origin of proteins involved (red and green for RefVirus and ProVirus, respectively). (D) The PHROGs Web site main page, where users can search for PHROGs or viruses of interest.

Figure 2 .

 2 Figure 2. Number of viruses considering their viral family and the class of their host (when not specified, viral families have a dsDNA genome). The size of the balloons are proportional to the number of viruses and the color reflects the proportion of proviruses.

Figure 3 .

 3 Figure 3. Cumulated number of clustered proteins. For example, point a means that for the standard clustering procedure, ∼234 000 proteins are in clusters that contain at least 200 proteins, whereas for the PHROG procedure, ∼390 000 proteins are in clusters >200 (point b). The inset at the top right is a zoom of the left part of the curve. The 5 largest PHROGs are highlighted by a cross at the bottom right (the two largest PHROGs at the bottom right gather 5795 and 5879 proteins).

Figure 4 .

 4 Figure 4. Identity percent (A) and coverage (B) for protein pairs in the same clusters. The clusters were separated according to interval of size, the first value 3 representing clusters that contain 3, 4, 5 or 6 proteins, 7 being clusters containing between 7 and 19 proteins, and the last interval 2000 being clusters >2000 proteins. Using the multiple alignments of each cluster, (i) the identity percent between two proteins is the number of amino acids that are identical in the two aligned proteins divided by the length of the smallest of the two proteins, and (ii) the coverage is the proportion of amino acids of one protein that is aligned to any amino acid (not to a gap) of the other protein. Subsamples of 1000 values where taken to draw each boxplot.

Figure 6 .

 6 Figure 6. Each vertex represents an annotation and two annotations are linked if considered as significantly colocalized (see Materials and Methods section). Edges were attributed a weight equal to their significance score. Only the 112 most frequent annotations (used for >650 proteins) are displayed.Annotations for which genes were colocalized with genes having the same annotation are drawn as squares. Among these 112, the 21 annotations not significantly colocalized to any other are displayed on the right.

Table 1 .

 1 Number of PHROGs (and proteins inside these PHROGs) annotated in the nine functional categories. The first line represent the whole dataset and the hierarchy in the following lines is represented by indentations. A color was attributed to each of the nine functional category

FUNDING

H2020 European Research Council [685778]. Funding for open access charge: Institut national de recherche pour l'agriculture, l'alimentation et l'environnement.

DATA AVAILABILITY

All results are accessible on the PHROG website available at https://phrogs.lmge.uca.fr/.

Sequences, multiple alignments, HMM profiles and annotations are downloadable as zipped archives on this web-site. An online manual helps users compare their protein datasets to the PHROG profiles.

SUPPLEMENTARY DATA

Supplementary Data are available at NARGAB Online. Conflict of interest statement. None declared.


Christophe Denis
email: christophe.denis@lip6.fr

Cadre méthodologique pour assurer la spécification d'un simulateur par apprentissage machine profond en vue de sa validation

Cadre méthodologique pour assurer la spécification d'un simulateur par apprentissage machine profond en vue de sa validation

Contexte et motivation industrielle

Le Groupe EDF exploite des matériels et des ouvrages à longue durée de vie comme des barrages hydroélectriques et des centrales nucléaires et mène historiquement au sein de sa R&D une forte activité de simulation numérique dans différents domaines scientifiques, allant de la mécanique des structures à la mécanique des fluides en passant l'électromagnétisme. Ces simulations ont pour objectif de faciliter la prise de décision industrielle sur des quesions de sûreté, de maintenance et d'optimisation de ses moyens de production électrique. A titre d'exemple, le code de simulation des écoulements à surface libre TELEMAC [ [START_REF] Hervouet | Hydrodynamics of Free Surface Flows: Modelling With the Finite Element Method[END_REF]], co-développée par la R&D d'EDF, est utilisé dans le but d'estimer la hauteur d'eau de fleuve afin de dimensionner les digues de protection des centrales nucléaires ou pour décider de l'emplacement d'une hydrolienne afin d'en optimiser sa production. Ce code de simulation comme une grande majorité d'autres codes industriels se base sur une modélisation hypothético-déductive, qui a permis d'obtenir des avancées significatives sur la prédiction et la compréhension de phénomènes. Le principe de cette modélisation est de formaliser le phénomène d'intérêt sous la forme d'un système d'équations mathématiques exprimées tout d'abord sur des supports continus puis discrétisées sous la forme d'un schéma numérique spatio-temporel. Une méthodologie éprouvée de vérification, de validation et de quantification d'incertitudes est mise en place à EDF R&D pour accréditer les résultats de ses codes de simulation.

La modélisation hypothético-déductive n'est pas toujours la panacée pour questionner un phénomène. En effet, il n'est pas toujours possible de trouver le modèle mathématique d'un phénomène mal connu dans l'état de connaissance scientifique actuel, d'une part, et la simulation de ce modèle peut produire des prédictions peu précises ou nécessiter des temps de calcul ou de consommation énergétique jugés trop importantes, d'autre part. Les techniques d'apprentissage machine profond rendent possible la mise au point de simulateurs sur des phénomènes pour lesquels on ne dispose pas de modèle théorique satisfaisant. Ainsi, de nombreuses disciplines scientifiques computationnelles s'orientent vers l'apprentissage machine profond pour mettre en oeuvre des simulateurs sans une formalisation préalable d'un modèle, par exemple pour la simulation de la turbulence en mécanique des fluides. Contrairement à l'approche hypothético-déductive, les simulateurs par apprentissage machine sont développés empiriquement sans être bordés par un cadre formel de spécification et de validation pénalisant leurs acceptabilités opérationnelles, et ce même si on peut vérifier le comportement empirique du simulateur en utilisant des techniques d'interprétation post hoc. Il peut être tentant faute de mieux d'adapter pour l'apprentissage machine profond la méthodologie éprouvée de spécification et de validation utilisée pour l'approche hypothético-déductive. Cette solution ne se justifie pas puisque la spécification et la mise au point d'une simulation par apprentissage machine sont différentes de celles mises en oeuvre pour une approche hypothético-déductive.

Nous proposons dans cette contribution un cadre méthodologique pour formaliser les spécifications d'un simulateur basée sur de l'apprentissage machine profond pour établir un cadre rigoureux de validation. Ce travail de recherche, mêlant épistémologie de la modélisation et informatique théorique, a pour objectif comme mentionné dans [START_REF] Girard-Satabin | Verification and validation of Machine Learning techniques[END_REF] de "réconcilier l'abondance des techniques de vérification de programmes classiques et l'absence de garanties sur les réseaux de neurones pour permettre aux logiciels critiques de conserver leur haut niveau de confiance". Notre contribution est organisée comme suit. Nous présentons tout d'abord les différences épistémologiques entre une modélisation hypothéticodéductive et une modélisation par apprentissage machine. Nous argumentons ensuite sur la nécessité de définir un cadre spécifique de validation de simulateurs par apprentissage, en adaptant à l'apprentissage machine la Théorie de la Modélisation et de la Simulationt pour border epistémologiquement l'utilisation de solutions logicielles de spécifications formelles. par rapport aux observations. Un modèle peut être une expérience de pensée comme celle de Galilée sur la chute des corps remettant en cause la loi aristotélicienne pour révéler le rôle joué par l'air. Galilée argumente également que l'on ne peut comprendre l'univers "si l'on ne s'applique d'abord à en comprendre la langue et à connaître les caractères avec lesquels il est écrit. • le nombre et la nature des données d'entrée X ainsi que la valeur de sortie y sont spécifiés par ces mêmes équations mathématiques.

On peut par exemple choisir comme modèle mathématique les équations de Saint-Venant pour étudier le comportement des écoulements à surface libre. L'évaluation numérique f (X) nécessite le développement d'un simulateur. Pour évaluer numériquement le comportement des écoulements à surface libre, le code de simulation Telemac, co-développé à EDF R&D, construit un schéma spatio-temporel discrétisé pour résoudre les équations mathématiques de Saint-Venant. L'acceptabilité opérationnelle du simulateur nécessite de confronter ces évaluations avec des observables du phénomène. Il est tout d'abord nécessaire de vérifier que le simulateur évalue précisément les trajectoires spécifiées par le modèle. Cette vérification est nécessaire mais non suffisante puisque l'on peut simuler des trajectoires très proches des spécifications du modèle et pour autant éloignées du comportement observable du phénomène. Sur le plan méthodologique, une démarche éprouvée de vérification, de validation et de quantification d'incertitudes a été proposée par [Obe10] :

• l'étape de vérification permet de contrôler la programmation du simulateur spécifiée par la fonction f ;

• l'étape de validation évalue la capacité de la fonction f à reproduire les observations du phénomène, ce qui revient à mettre à l'épreuve les spécifications de la fonction f ;

• l'étape de quantification des incertitudes permet d'estimer des intervalles de confiance par rapport aux évaluations de la fonction f , pour tenir compte d'un manque de connaissance de certains comportements du phénomène et des incertitudes sur les paramètres d'entrée.

Cette démarche de vérification et de validation est basée sur l'hypothèse d'une séparation des causes à l'origine de l'écart entre les prédictions et les observations. Cette hypothèse est en cohérence avec les processus de conception et de développement d'un simulateur basée sur une modélisation hypothéticodéductive étant donné que le développement du simulateur est toujours encadré par le choix préalable d'un modèle mathématique, qui peut être à son tour remis en cause durant l'étape de validation.

La modélisation hypothético-déductive n'est pas toujours la panacée pour questionner un phénomène. En effet, il n'est pas toujours possible de trouver un modèle mathématique suffisamment précis pour décrire un phénomène mal connu dans l'état de connaissance scientifique actuel, d'une part, et la simulation de ce modèle peut produire des prédictions peu précises ou nécessiter des temps de calcul ou de consommation énergétique jugés trop importantes, d'autre part. Ainsi, de nombreuses disciplines scientifiques computationnelles s'orientent vers l'apprentissage machine profond pour mettre en oeuvre des simulateurs sans une formalisation préalable d'un modèle.

Simulation par apprentissage machine profond

Les techniques d'apprentissage machine profond rendent possible la mise au point de simulateurs sur des phénomènes pour lesquels on ne dispose pas de modèle satisfaisant. En reprenant la représentation fonctionnelle de la figure 1, le simulateur prédit une valeur y en fonction de paramètres d'entrée X.

Conclusion et perspectives

  Il est écrit dans la langue mathématique et ses caractères sont des triangles, des cercles et autres figures géométriques, sans le moyen desquels il est humainement impossible d'en comprendre un mot. Sans eux, c'est une errance vaine dans un labyrinthe obscur". Initiée par Galilée et Descartes, la modélisation hypothético-déductive a permis d'obtenir des avancées significatives sur la prédiction et la compréhension de phénomènes, non nécessairement physiques d'ailleurs, ce que Wigner qualifiera d'efficacité déraisonnable. Le principe de cette modélisation est de prendre de la distance par rapport aux observations obtenues sur le monde sensible, en formalisant le phénomène sous la forme d'un système d'équations mathématiques qui sont le plus souvent exprimées sur des supports continus. Nous représentons le modèle par une fonction f , comme présenté à la figure 1.

Figure 1 :

 1 Figure 1: Représentation fonctionnelle d'un modèle Dans le cadre de l'approche hypothético-déductive,

Figure 2 :

 2 Figure 2: Représentation fonctionnelle d'un simulateur par apprentissage machine La mise au point du simulateur ne découle pas dans la majorité des cas d'une spécification de la fonction f et des données d'entrée comme dans le cadre de la modélisation hypothético-déductive. Cette spécification s'établit implicitement lors de la mise au point du simulateur. La conception du simulateur suit une

Figure 3 :

 3 Figure 3: Synthèse visuelle de la méthodologie empirique de mise au point d'un simulateur par apprentissage machine Ainsi, au fil des itérations i de cette boucle de rétroaction, une succession de fonctions f i et de représentations des données d'entrée X i est produit jusqu'à ce qu'un seuil sur la métrique d'erreur soit obtenue avec une méthode d'optimisation. Il est montré dans [NPS21] que les techniques d'optimisation des réseaux de neurones profonds introduit une rupture avec le phénomène cible, pouvant conduire à des spécifications de problèmes différentes de celui du cahier des charges. Ce glissement de spécification s'illustre dans l'exemple bien connu d'un projet de reconnaissance automatique d'un loup ou d'un husky dans des photographies. L'utilisation d'une méthode d'interprétabilité post-hoc a montré que la spécification du simulateur obtenu était différente du cahier des charges car détectant ou non la présence de neige sur les photos. Pour se prémunir en particulier de ce glissement, nous souhaitons mettre en place un cadre formel de spécification pour appliquer systématiquement une méthodologie de vérification et de validation, comme c'est la cas pour les simulateurs basés sur une approche hypothético-déductive.

Table 1 :

 1 Cette méthodologie consiste à définir les spécifications formelles du simulateur par apprentissage machine, son modèle sous-jacent, afin de mettre en place ensuite en place une démarche rigoureuse de validation. Le terme modèle est par nature polysémique et peut prêter à confusion puisque que ce terme n'a pas la même signification en science des données et en épistémologie de la simulation. Il n'existe de distinction en science des données entre le simulateur (dénommé confusément par le terme modèle) qui s'exécute sur un ordinateur, et ses spécifications. Le terme modèle est considéré ici comme un objet médiateur entre un observateur et un phénomène tel que caractérisé par Marvin Minsky "Pour un observateur B, un objet A* est un modèle d'un objet A, dans la mesure où B peut utiliser A* pour répondre à des questions qui l'intéressent au sujet de A"[START_REF] Minsky | Matter, mind, and models[END_REF]. Le tableau 1 présente une synthèse du cadre méthodologique de spécification et de validation pour un simulateur construit selon une approche de modélisation hypothético-déductive ou par apprentissage machine profond. Synthèse du cadre méthodologique de spécification.La caractérisation de Minsky ne peut pas être utilisé pour formaliser les spécifications car elle distingue pas le simulateur du modèle et ne prend pas en compte la méthodologie empirique de mise au point du simulateur. Par contre, il est possible de définir un cadre épsitémologique pertinent en utilisant les concepts issus de la Théorie de la Modélisation et de la Simulation (TMS) proposée dans les années 1970 par Bernard P. Ziegler[START_REF] Bernard | Theory of Modeling and Simulation-Integrating Discrete Event and Continuous Complex Dynamic Systems[END_REF].

	4 Cadre méthodologique pour la spécification de
	simulateurs basés sur de l'apprentissage ma-
	chine profond
	Le cadre méthodologique que nous proposons s'inspire des travaux menés par
	[Que12] pour spécifier des simulations multi-agent, d'une part, et se base sur la
	clarification épistémologique de l'apprentissage machine proposée par [DV22],

partir de son cahier des charges, d'une part, et de quantifier l'apport de mécanismes d'interprétabilité post-hoc afin de le comparer avec le niveau souhaité par le cahier des charges. Il est enfin nécessaire d'associer à ce cadre un travail de fond en informatique théorique pour proposer des nouveaux modèles de calcul pour l'apprentissage machine dépassant la thèse de Church-Turing

[START_REF] Shagir | The Nature of Physical Computation[END_REF]

.

  Les techniques d'apprentissage machine profond permettent de concevoir des simulateurs performants sur des phénomènes ne disposant pas d'un modèle théorique satisfaisant. Cette absence de spécification formelle empêche de mettre en place une démarche rigoureuse de validation. En effet, les paramètres d'un simulateur par apprentissage machine profond s'ajustent dans une boucle de rétroacation visant à réduire une métrique d'erreur. L'optimisation des paramètres est agnostique par rapport au phénomène d'intérêt pouvant induire un glissement implicite de spécification par rapport à celle souhaitée dans le cahier des charges. Ce glissement de spécification peut être détecté en appliquant des techniques d'interprétabilité post-hoc lors d'analyses de dysfonctionnements du simulateur. Cette stratégie n'est pas acceptable pour des applications industrielles suffisamment critiques. Cette contribution constitue le fondement épistémologique d'une démarche systématique de validation ne reposant pas sur l'hypothèse peu crédible et non utile d'une compréhension fine des mécanismes internes du simulateur. Il s'agit d'expliciter les spécifications formelles du simulateur, son modèle sous-jacent, en adaptant pour l'apprentissage machine profond la théorie de la Modélisation et de la Simulation, permettant de border l'utilisation d'une solution logicielle de spécification formelle de simulateurs basées sur de l'apprentissage machine profond. Ce cadre méthodologie est en cours de mise à l'épreuve sur un cas d'usage prédisant un phénomène physique en mécanique des fluides et sur la production d'un résumé de texte avec l'agent conversationnel ChatGPT. Enfin, nous poursuivons un travail de clarification épistémologique pour définir la notion d'algorithme d'apprentissage machine associé au simulateur en nous basant sur des nouveaux modèles de calcul en informatique théorique.

Modélisation hypothético-déductiveL'acquisition de la connaissance sur un phénomène nécessite le plus souvent la mise en place d'un objet médiateur que l'on appelle modèle, pour prendre recul