Etude des séquences modèles de RdRps

Nous avons constitué un jeu de données modèles de RdRps qui nous servira pour les jeux d’apprentissage. Pour ce jeu de données, les polymérases ont été obtenues sur Uniprot en ne gardant que les séquences virales de Swissprot ayant un domaine RdRp. Cette sélection cor-respond à 599 séquences. Puis, les séquences ont été filtrées en fonction de la taxonomie dans le but d’éliminer les groupes de virus qui infectent uniquement des mammifères par exemple (comme le groupe des virus de la grippe). En effet, les virus d’intérêt pour cette étude, sont des virus infectants plutôt les micro-organismes du plancton (micro-algues, bactéries, etc.). Après ce premier filtrage, réalisé à l’aide du ICVT [11], le nombre de séquences est des-cendu à 410. La liste des séquences sélectionnées en fonction de leur taxonomie est acces-sible sur http://www.irisa.fr/dyliss/public/Virus/mjusot/selection_ taxonomie.xls. Etant donné que ces séquences modèles vont être la base de cette recherche, il est important de ne garder uniquement que les données les plus fiables, préférentiellement dis-posant de preuves expérimentales. Ainsi, un second filtrage a été réalisé en fonction du score d’annotation Uniprot. Celui-ci repose sur une mesure heuristique de la qualité d’annotation et correspond à un chiffre entre 0 et 5 (5 étant le score maximal). Seules les séquences possédant au minimum un score d’annotation de 4 ont été gardées pour la suite soit un total de 288 séquences. Enfin, la redondance des séquences a été réduite à 70% par CD-hit [15] ce qui nous a permis d’obtenir au final 124 séquences de RdRps de divers virus à ARN.

Un nettoyage a également été effectué sur ces données en enlevant les séquences conte-nant des acides aminés inconnus ’X’. En effet, Protomata-Learner n’accepte pas les séquences contenant ce type d’acides aminés pour apprendre les modèles. Cette étape a été inclue dans l’automatisation de la procédure. Nous avons finalement obtenu 108 séquences pour notre jeu de données modèles.

Pour tester les modèles, il est nécessaires d’avoir des jeux de validation. Pour cela, différents jeux de validation ont été constitués :

– un premier utilisé pour la validation croisée contenant 27 séquences de polymérases sélec-tionnées sur SwissProt et ayant un score d’annotation d’au moins 4 mais n’étant pas des RdRps virales (6 RdRps eucaryotes, 13 polymérases d’archées et 8 polymérases de virus à ADN). Ces séquences dites « négatives » mais relativement proches de ce que l’on recherche ont été sélectionnées pour aider à trouver les paramètres les plus discriminants entre les sé-quences positives (RdRps virales) et celles-ci. On rajoutera en effet, à chaque étape du k-fold, m séquences positives du jeu modèle de RdRps pour former le jeu de validation.

– Un second jeu de données de validation a été obtenu en prenant toutes les séquences de Swissprot et en ajoutant également les séquences Uniprot étant annotées RdRp. Ce jeu a été nettoyé pour retirer toutes les séquences contenues dans le jeu d’apprentissage. Il a également été filtré par Uniref50. Ce jeu de validation nous servira pour une comparaison des résultats entre les 3 méthodes de caractérisation des RdRps.

Pour générer une banque de capsides virales, toutes les séquences de la Protein Data Bank (PDB) contenant le terme Gene Ontology « Viral Capsid » ont été récupérées, ainsi que toutes les structures annotées capsides dans la Viral Protein Structural Database [16]. Les deux listes de structures ont été fusionnées puis nettoyées de sorte à enlever toutes les chaînes qui n’étaient pas des protéines. Une banque non redondante a alors été créée en utilisant Uclust (version 8.0.1623) [17]. A partir de ces résultats, pour chaque cluster, la protéine de meilleure résolution et préférentiellement résolue par cristallographie a été choisie comme représentant. La banque finale contient un total de 327  chaînes PDB.

Un fragment de contact (CF) est défini comme une paire de segments de squelette carboné d’une protéine qui sont proches dans l’espace (structure 3D) ou plus précisément qui partagent un contact (distance entre C inférieure à =7 Angstroms), et sont suffisamment proches pour interagir (chaque C est au plus loin à =13 Angstroms d’un C de l’autre segment) (voir figure 3) [12]. Pour extraire ces CFs, nous avons utilisé un script, fourni par C. Galiez, auquel nous avons fourni la liste des chaînes PDB correspondant aux structures de capsides précédemment sélectionnées. Ce script nous permet ainsi d’obtenir pour chaque chaîne les CFs lui correspondant sous format PDB.

Beaucoup de CFs sont générés et ils ne sont pas forcément spécifiques des capsides tant en structure qu’en séquence. Nous avons voulu dans un premier temps essayer de déterminer quels sont les CFs ou sous ensembles de CFs qui sont spécifiques des capsides virales. Pour cela, nous avons voulu rechercher les CFs extraits de chaque capside dans la PDB. Notre choix s’est porté sur Yakusa [13] qui est un programme permettant de rechercher dans une banque toutes structures protéiques similaires à une structure requête car la recherche est rapide, locale et sans gap ce qui est particulièrement adapté à notre cas. De plus comme le code est disponible, il a été facile de l’adapter à nos besoins.

Yakusa transforme toutes les structures (requête et protéines de la banque) en enchaînement d’angles internes (angle soit l’angle dihédral entre 4 carbones alpha consécutifs). La similarité entre séquences d’angles est établie selon les étapes suivantes :

1. Construction d’un automate déterministe fini décrivant tous les mots structuraux d’une taille donnée (souvent 4) de la protéine requête, c’est à dire toutes les combinaisons de taille fixée d’angles alpha existants.

2. Recherche de tous les mots identiques ou proches de ceux de la protéine requête dans chaque structure de la banque.

FIGURE 4 – Définition des filtres sur les longueurs et sur les distances ajoutés dans Yakusa. A gauche en rouge, CF requête et à droite en bleu, couple de SHSPs trouvé par Yakusa corres-pondant au CF requête. La partie colorée en bleue du CF est la partie correspondant aux SHSPs trouvés. Les longueurs L1a et L2a correspondent aux longueurs des 2 segments du CF requête, tandis que les longueurs L1b et L2b correspondent aux longueurs des 2 SHSPs associés. Pour être considérés comme compatible la longueur L1b (respectivement L2b) doit être supérieure ou égale à n % de L1a (respectivement L2a), n étant fixé par défaut à 80%. Les distances dia, dma et df a sont comparées avec les distances dib, dmb et df b. La distance dib doit entre égale à dia à plus ou moins m % (m étant fixé par défaut à 10%). Il en est de même entre dmb et dma ainsi qu’entre df b et df a.

3. Filtres des mots répétés et extension des mots trouvés vers des sous-structures (SHSPs pour structural high-scoring pairs) plus longues.

4. Sélection de SHSPs compatibles pour chaque protéine de la banque.

5. Classement des structures basé sur 3 scores : similarité, compatibilité et probabilité.

Ce programme a été adapté ici à partir du code existant de Yakusa en C. La banque utilisée est la PDB non redondante 50% du site de la PDB datant de mai 2015. Pour chaque capside, un fichier PDB contenant tous les CFs identifiés a été généré. Chaque CF est considéré comme une chaîne différente qui est en deux fragments. Chacun de ces fichiers PDB de CFs est utilisé comme requête. On obtient alors un ensemble de SHSPs correspondant chacun à un fragment de CFs. Les deux dernières étapes de Yakusa ont été adaptées afin de filtrer les SHSPs trouvés non compatibles avec la définition des CFs. Nous avons donc mis en place deux filtres qui remplacent la sélection des SHSPs compatibles :

– Le premier filtre repose sur la longueur du SHSP trouvé par rapport à la longueur du segment de CF requête. Le SHSP trouvé doit être de longueur supérieure ou égale à un pourcentage de la longueur de la requête. Ce pourcentage seuil a été fixé par défaut à 80% mais a été rajouté dans les paramètres de Yakusa est peut ainsi être modifié éventuellement dans l’optique de trouver les sous-fragments conservés qui pourraient être plus caractéristiques des structures.

– Le second filtre repose sur la recherche de paires de SHSPs qui doivent être compatibles en distances pour être considérés comme des CFs. Pour cela, il faut déjà s’assurer que les 2 SHSPs proviennent chacun d’un fragment différent du même CF de la protéine requête puis calculer les distances entre les SHSPs et s’assurer qu’elles sont à un pourcentage près les mêmes que celles du CF requête. Ce seuil a été fixé à 10% par défaut mais a également été rajouté en paramètre. Pour être rapide, les distances testées correspondent au début du SHSP, au milieu et à sa fin (voir figure 4).

Nous disposons des données métagénomiques de Tara Oceans des stations 122 à 125 (sta-tions des îles Marquises). Les données utilisées sont un assemblage des gènes (issus d’un sé-quençage NGS) correspondant au mélange des métatranscriptomes de toutes les stations (soit 53 959 998 séquences toutes fractions confondues entre 0,8  m et 2000 m).

Après toutes les étapes de nettoyage, le jeu final d’apprentissage pour les RdRps contient 108 séquences réparties de la manière suivante dans la taxonomie des virus : 2 séquences de virus à ARN double brins, 12 séquences de virus à ARN simple brin à polarité négative et 94 séquences de virus à ARN simple brin à polarité positive. La longueur de ces séquences est très variable (entre 665 et 3945 résidus) avec une longueur moyenne de 2340 acides ami-nés par séquences et un écart type de 704,6 (distribution présentée figure figA1 en annexe). Cette distribution n’est pas surprenante car la majeure partie de ces séquences virales sont des polyprotéines, c’est à dire le produit d’expression d’un seul gène en une unique protéine géné-ralement non fonctionnelle, permettant ainsi aux virus d’avoir un génome très compact. Cette unique protéine est ensuite clivée en plusieurs protéines plus petites qui sont alors fonction-nelles. Ce point est important à souligner car le fait de caractériser ces polyprotéines, et non uniquement la partie de séquences correspondant aux RdRps, peut biaiser l’analyse. En effet, une séquence d’hélicase peut être reconnue par BLAST par exemple parce qu’elle sera proche de la partie hélicase d’une polyprotéine. Toutefois, si la conservation est très faible dans les parties non RdRps, les méthodes vont avoir tendance à se focaliser sur la partie plus conservée. Notamment en ce qui concerne Protomata, le programme devrait s’affranchir de cette barrière en se concentrant en particulier sur les blocs conservés sur l’ensemble ou sur une majorité de séquences et non sur les parties conservées dans les sous groupes de virus.

Des tests effectués préalablement à la validation croisée ont été réalisés sur des sous en-sembles de données (petits jeux issus des alignements PFams des différentes familles de RdRps). Ces tests ont permis non seulement de guider le choix des paramètres à tester vers certaines valeurs (et ainsi éviter d’en tester trop) mais également de confirmer le choix de certains para-mètres fixés comme par exemple la mixture de Dirichlet (voir Figure A3).

Pour la validation croisée, le k-fold choisi a été de 3. Prendre un k plus petit aurait rendu la procédure plus longue. De plus, ceci permet d’ajouter m=36 séquences positives au jeu de validation, ce qui est relativement d’un même ordre de grandeur que le nombre de séquences du jeu négatif, et permet ainsi de ne pas trop biaiser les résultats statistiques de la validation croisée. Les différentes valeurs des paramètres testés sont présentées dans la table 1.

Paramètres Valeurs testées

Un nombre total de 360 combinaisons de paramètres ont finalement été testées. Ces expériences ont été parallélisées par groupe de 12 à 16 combinaisons pour optimiser le temps de calcul. Pour chaque combinaison de paramètres, les valeurs suivantes ont été calculées en fonctions de 2 seuils différents pour la prédiction (voir ci dessous) :

– la distribution des scores pour chaque k-fold, en séparant les séquences positives des néga-tives, permettant ainsi de voir quels paramètres les discriminent le mieux pour l’ensemble des k-folds.

– la marge correspondant à la différence entre la moyenne des scores des séquences positives et la moyenne des scores des séquences négatives ramenée à la somme de leur écarts-types. Ceci permet également de voir quels paramètres les discriminent le mieux pour l’ensemble des k-folds. M arge = (moyenne(scores_positifs)  moyenne(scores_negatifs)) ecart  type(scores_positifs)+ecart  type(scores_negatifs)

– le rappel mesurant le nombre de prédictions positives (V rai_P ositif s) parmi les séquences de RdRps (V rai_P ositif s + F aux_N egatif s). Rappel = V rai_P ositifs (V rai_P ositifs+F aux_Negatifs)

– la précision mesurant le nombre de séquences de RdRps (V rai_P ositif s) parmi les sé-quences prédites positives (V rai_P ositif s + F aux_P ositif s). P recision = V rai_P ositifs (V rai_P ositifs+F aux_P ositifs)

Les trois méthodes de caractérisation des RdRps ont été appliquées avec pour jeu d’ap-prentissage, le jeu complet de séquences modèles et pour jeu de validation le jeu basé sur les séquences de Swissprot. Pour chaque méthode, la courbe précision-rappel a été calculée ainsi que les valeurs des tables de contingences pour différents seuils. Il faut noter que le jeu de validation contient très peu de RdRps. L’utilisation de la courbe précision-rappel est adapté à ce genre de données présentant un biais [18] (contrairement aux analyses de courbes ROC par exemple). BLAST Cette méthode a été utilisée avec un seuil de e-value assez lâche fixé à 0,1. L’aire sous la courbe précision-rappel (figure 7 (a)) est à 0,72. On peut voir 3 pentes sur la courbe, ce qui laisse à penser qu’il y a 3 coupures concernant les e-values. La première coupure est immédiate : on passe d’une précision de 1 à 0,96 dès les premières valeurs car des faux posi-tifs (séquences prédites RdRps mais n’en étant pas) sont immédiatement reconnus par BLAST comme étant des RdRps. En effet, 6 faux positifs sont immédiatement détectés. Cinq de ces sé-quences correspondent à des polyprotéines de virus à ARN simple brin n’étant pas des RdRps mais contenant des protéines de capside, d’enveloppe, des protéases ou encore des hélicases. Le dernier faux positif correspond à une RdRp eucaryote. La deuxième pente de la courbe (chute de la précision pour un rappel supérieur à 0,75) correspond aux séquences qui n’ont pas été reconnues par BLAST comme étant des RdRps virales. Le logiciel ne leur a ainsi pas attribué de e-value et pour pouvoir réaliser cette courbe nous leur en avons mis une par défaut à 1 (donc supérieur au seuil), ce qui explique cette chute. Et enfin, la troisième petite pente correspond à un artéfact de la méthode de calcul de la courbe précision rappel. Celle-ci va en effet récupérer à nouveau des vrais positifs parmi les séquences qui n’ont pas été reconnues ce qui créé une petite augmentation locale de la précision (pour un rappel entre 0,80 et 0,90) mais qui ne doit pas être prise en compte.
Les tables de contingence sont présentées table 2. On constate que si l’on prend une e-value plus stringente à 10 4, on perd de nombreux vrais positifs (441) pour une perte minimale de faux positifs (9). En effet, le rappel passe de 0,74 à 0,50 et même la précision passe de 0,88 à 0,83. Ainsi, ce choix ne semble pas plus approprié. Parmi les faux positifs trouvés avec une e-value à 0,1, il y a 127 séquences eucaryotes dont 44 correspondant à des RdRps eucaryotes, ainsi que diverses enzymes. On retrouve également 58 protéines de virus correspondant essentiellement à des hélicases et protéinases.
HMMER Cette methode a également été utilisée avec un seuil de e-value assez lâche fixé a 0,1. L’aire sous la courbe précision-rappel (figure 7 (b)) est à 0,44. On peut voir comme pour BLAST les 3 mêmes pentes sur la courbe. La première pente descend plus vite que BLAST mais arrive un peu plus tard. En effet, on obtient le premier faux positif pour une e-value de 10 194 alors que l’on a déjà trouvé 70 vrais positifs. La fin de cette première pente se situe autour d’une e-value de 10 98. d’une précision de 0,98 à 0,89. Ces 19 faux positifs sont presque tous des protéinases de virus à ARN simple brin à l’exception de 3 polyprotéines de virus à ARN simple brin contenant des séquences de protéines d’enveloppe et de capside ainsi que des protéases. La deuxième pente arrive pour un rappel légèrement supérieur à 0,40 soit beaucoup plus tôt que pour BLAST. Et on retrouve le même artéfact qu’avec blast. Les tables de contingence sont présentées table 3. De même que pour BLAST, on constate que si l’on prend une e-value plus stringente à 10 4, on perds 102 vrais positifs pour seulement 7 faux positifs. On passe alors d’un rappel à 0,43 à une valeur de 0,38 pour une même précision à 0,94. Ainsi, ce seuil plus stringent ne semble également pas adapté. Concernant les 51 faux positifs trouvés avec une e-value à 0,1, 20 séquences correspondent à des RdRps eucaryotes, mais on a également plusieurs enzymes ou protéines eucaryotes présentant un site de fixation à l’ARN ou l’ADN.
Protomata Pour cette méthode, il a fallu tester différents paramètres p pour trouver celui adapté au nombre de séquences. Ainsi, nous avons gardé ceux présentant les meilleures aires sous la courbe précision-rappel (voir figure A4 et A5 en annexe). Ainsi pour le jeu de paramètre 1.2.1, nous avons choisi p = 10 ayant une aire sous la courbe précision-rappel à 0,65 (figure 7 (c)). Tandis que pour le jeu 8.1.4, nous avons choisi p = 5 qui présente une aire sous la courbe à 0,68 (figure 7 (d)). Pour les deux courbes précision-rappel, on retrouve 2 pentes : la première ressemblant à celle trouvée chez BLAST et HMMER correspond à l’identification de 14 faux positifs ayant de très bons scores et correspondant tous à des protéinases de virus à ARN simple brins. Les deux modèles trouvent les même faux positifs lors de cette pente. La deuxième pente correspond à une augmentation des faux positifs. Elle commence plus tôt pour le jeu 1.2.1 (à savoir pour un plus faible rappel à 0,45 environ contre 0,6 pour le jeu 8.1.4). Il n’y a pas d’artéfact sur ces courbes puisque Protomata a calculé un score pour chaque séquence. Les courbes précision-rappel nous permettent ici de fixer un seuil de score plus adapté que ceux utilisés pour la validation croisée. En effet, il faut trouver les seuils de scores qui optimisent aussi bien le rappel que la précision et pour cela, on peut choisir de se placer sur la courbe juste avant la deuxième pente et trouver les seuils de scores associés à ces valeurs. Ainsi, pour jeu de paramètre 1.2.1, on a choisi un seuil à 21 permettant d’obtenir un rappel à 0.53 et une précision à 0,87. Et pour le jeu 8.1.4, on a choisi un seuil à 35 permettant d’obtenir un rappel à 0.58 et une précision à 0,93. Les tableaux de contingence associés sont présentés table 4.
Dans les 2 cas, les faux positifs trouvés sont essentiellement des protéinases et hélicases de virus à ARN, des RdRps eucaryotes, des enzymes eucaryotes présentant un site de fixation à l’ARN ou l’ADN (par exemple la Ribonuclease J). Le modèle 1.2.1 a également trouvé une reverse transcriptase de retrovirus (protéine P11283).
Discussion En comparant les valeurs statistiques obtenues avec les différentes méthodes, on constate que BLAST présente de meilleurs résultats, il trouve en effet beaucoup plus de vrais positifs sans augmenter énormément le nombre de faux positifs. Ces résultats sont surprenants puisqu’il était attendu que HMMER et Protomata soient plus appropriés pour les séquences vi-rales qui présentent un fort taux de mutation. Concernant HMMER, il est possible au vu de sa courbe précision-rappel 7(b), que la e-value soit plus stringente que BLAST. Il faudrait envi-sager de regarder des seuils encore plus lâches que 0,1. Concernant Protomata, malgré tout le travail de paramétrisation effectué, cet outil ne parvient pas à faire mieux que BLAST sur les polymérases. Ainsi, cet outil n’est peut être pas adapté à l’étude des séquences virales diver-gentes. Concernant les faux positifs, les 3 méthodes retrouvent le même type de séquences : des séquences de virus à ARN non RdRps, des RdRps eucaryotes et d’autres séquences eucaryotes. Ces résultats sont à la fois surprenants et attendus. En effet, dans la littérature les RdRps virales et eucaryotes sont décrites comme partageant peu d’homologie de séquences [19] [8]. Il est donc surprenant de les retrouver aussi facilement. Toutefois, le motif DxDGD, essentiel pour l’activité des RdRps cellulaires, semble être un reste du motif GDD des RdRps virales [19]. Cela met donc en avant une certaine conservation entre les RdRps eucaryotes et celles virales, même si elles appartiennent à des super-familles différentes et cela confirme la difficulté à les distinguer. Egalement, on trouve selon les méthodes, des séquences plus ou moins nombreuses d’enzymes eucaryotes présentant un site de fixation à l’ADN ou l’ARN. Ainsi peut être que la similarité de séquence se situe essentiellement aux niveaux des domaines de fixation à l’ARN ou l’ADN. A l’inverse, un résultat prévisible est le fait de trouver des séquences de virus à ARN n’étant pas des RdRps. Ceci est très certainement dû au fait de travailler sur des poly-protéines comme expliqué dans la partie sur l’étude des séquences modèles de RdRps. Ainsi, même Protomata n’est pas capable de passer cette barrière et pour améliorer les résultats il serait nécessaire de cliver les séquences du jeu modèle pour ne garder que la partie correspondant aux RdRps. Cela permettrait assurément de diminuer le nombre de faux positifs en éliminant toutes les protéinases, hélicases et autres protéines virales identifiées. Et cela permettrait ainsi de faire disparaître la première pente de la courbe précision-rappel. Il est possible que cela influe éga-lement sur la distribution des scores de Protomata et améliore ainsi la discrimination entre les séquences négatives et positives.

Après extraction des CFs sur l’ensemble de la banque de structures de capsides, 11395 CFs ont été obtenus. Les chaînes de structures utilisées présentent entre 1 et 168 CFs avec une moyenne de 25 CFs par chaîne et un écart-type de 26,9 (voir distribution en annexe figure A6). Les fragments de CFs présentent une longueur moyenne de 16,9 résidus avec un écart-type de 11 (voir figure A7 en annexe). Il faut comprendre que 2 fragments formant un CF n’ont pas forcément la même longueur et également que ceux-ci peuvent être chevauchants. Ainsi, certains CFs sont très longs (une dizaine de fragments de CFs dépassent les 100 résidus). Ce sont des fragments qui s’alignent avec eux même (auto-fragments). Mais ces cas sont rares puisque 89% des fragments de CFs sont d’une longueur inférieure à 30 résidus. Dans 12 chaînes de capsides, aucun CFs n’a été trouvé.