Un promoteur U6 (mU6) de souris élargit les sites de ciblage génomique
Une conception optimale de bibliothèque d’ARNg nécessite un grand nombre de sites cibles génomiques accessibles. La transcription de l’ARNg est généralement pilotée par le promoteur humain U6 (hU6), dont on pense qu’il nécessite la guanine (G) comme premier nucléotide de sa transcription1,2,3. Dans le cas où le premier nucléotide n’est pas un G, il est possible de remplacer le premier nucléotide par un G ou d’ajouter un G supplémentaire à l’extrémité 5′ de l’ARNg, ce qui entraîne un mismatch ARNg-ADN à l’extrémité 5′. WT-SpCas9 peut tolérer les mésappariements ARNg-ADN à l’extrémité 5′, il peut donc cibler toute séquence N20NGG avec le promoteur hU61,2,3. Cependant, les nucléases Cas9 hautement spécifiques telles que eSpCas9(1.1) et SpCas9-HF1 sont sensibles aux mésappariements ARNg-ADN à l’extrémité 5′29. Elles ne peuvent cibler la séquence GN19NGG que lorsque le promoteur hU6 est utilisé, ce qui limite la sélection du site cible.
Une étude précédente a montré que le promoteur U6 (mU6) de la souris peut initier une transcription en adénine (A) ou en G30, ce qui pourrait potentiellement élargir la sélection des cibles. Nous avons comparé l’activité du promoteur mU6 et du promoteur hU6 par l’expression transitoire de gRNA pour l’édition du génome dans des cellules HEK293T exprimant WT-SpCas9 (Fig. 1a). Le promoteur mU6 a montré une activité similaire à celle du promoteur hU6 pour 12 gRNA testés initiés par G (Fig. 1b, c, Tableau 1 ; Fig. 1a supplémentaire). Nous avons testé neuf ARNg initiés avec A mais sans G sur 1-4 nucléotides, évitant ainsi les ARNg tronqués fonctionnels transcrits à partir de G (Tableau 1). À notre grande surprise, les deux promoteurs pouvaient favoriser l’édition du génome avec des ARNg initiés par A (Fig. 1d, e ; Fig. 1b supplémentaire). Parmi les neuf ARNg testés, seul un ARNg (A8) piloté par hU6 a montré une faible efficacité. Nous avons testé neuf autres ARNg initiés par A, mais contenant G à 1-5 nucléotides. Dans ce cas, le promoteur mU6 a montré une activité généralement plus élevée que le promoteur hU6 (figures supplémentaires 2a-c). Ensuite, nous avons comparé l’activité du promoteur mU6 et du promoteur hU6 dans les cellules HeLa, et ils ont montré une activité similaire (figures supplémentaires 3a-f). De plus, nous avons comparé l’activité du promoteur mU6 et du promoteur hU6 dans un vecteur lentivirus, et ils ont montré une activité similaire à deux moments, le troisième et le cinquième jour, après la transduction (figures supplémentaires 4a, b). Nos résultats sont cohérents avec une étude très récente selon laquelle le promoteur hU6 peut transcrire des petits ARN initiés avec A31.
Les deux promoteurs mU6 et hU6 permettent de transcrire des gRNA initiés par A ou G pour l’édition du génome. a Schéma de comparaison entre mU6 et hU6. Ces deux promoteurs ont été utilisés pour transcrire des ARNg initiés avec un nucléotide A ou G. b, c Comparaison des promoteurs mU6 et hU6 pour l’édition du génome avec des ARNg initiés avec G. Les données sont présentées sous forme de moyenne ± s.d. (n = 2). d, e Comparaison du promoteur mU6 pour l’édition du génome avec des ARNg initiés avec A. Les données sont présentées sous forme de moyenne ± écart-type (n = 2). f Comparaison du promoteur mU6 pour l’édition du génome avec des ARNg initiés avec A ou G trois jours après la transfection. g Comparaison du promoteur mU6 pour l’édition du génome avec des ARNg initiés avec A ou G cinq jours après la transfection. Les données sont présentées sous forme de moyenne ± écart-type. P > 0,05 ; P < 0,05 par ANOVA à deux voies (n = 2). Les données sources sont fournies sous forme de fichier Source Data
Puis, nous avons testé l’activité du promoteur mU6 pour l’édition du génome avec des gRNA initiés par C ou T, mais modifiés en A ou G (Tableau 1). Après 3 jours d’édition du génome, les ARNg initiés avec G ont montré une activité plus élevée pour trois des cinq ARNg testés (Fig. 1f ; Fig. 5a supplémentaire), mais la différence a été éliminée après 5 jours (Fig. 1g ; Fig. 5b supplémentaire). Le promoteur mU6 a été choisi dans l’étude suivante.
Une stratégie de test à haut débit de l’activité des ARNg
Une étude récente a montré qu’une stratégie de paire ARN guide-cible permet de tester à haut débit l’activité des ARNg pour Cpf132. Dans cette stratégie, les séquences guide-ARN-cible synthétisées sont délivrées dans des cellules exprimant Cas9 par des lentivirus (Fig. 2a). Après édition du génome, les séquences cibles sont amplifiées par PCR pour un séquençage profond, ce qui permet de mesurer directement les taux d’insertion/délétion (indel) induits par les nucléases Cas9. Un autre avantage est que les lentivirus s’intègrent préférentiellement dans les régions transcriptionnellement actives qui sont beaucoup plus accessibles pour la machinerie CRISPR/Cas932,33,34, minimisant l’influence de l’édition du génome par l’accessibilité de la chromatine. Par conséquent, l’ensemble des données obtenues par cette stratégie offre la possibilité d’élucider l’activité inhérente des ARNg en se basant exclusivement sur leurs caractéristiques de séquence.
Une stratégie de paire ARN- cible guide pour le criblage à haut débit de l’activité des ARNg dans les cellules humaines. a Schéma de la stratégie de paire ARN- cible guide pour le test d’activité des ARNg. Un vecteur lentiviral contient un promoteur mU6 et une paire guide-ARN-cible. Le vecteur a été utilisé pour transduire des cellules exprimant les nucléases Cas9. Les indels seraient induits au niveau des cibles intégrées par les ARNg correspondants. b Schéma de la conception et du criblage à haut débit d’une bibliothèque d’ARNg. Une bibliothèque de 80 263 paires ARN guide-cible a été conçue et synthétisée par microarray. Les oligonucléotides ont été amplifiés par PCR et clonés dans des vecteurs lentiviraux par assemblage de Gibson. La bibliothèque a été emballée dans des virus et transduite dans des cellules exprimant les nucléases Cas9 pour l’édition du génome. Les sites cibles intégrés ont été amplifiés par PCR pour une analyse de séquençage profond. c La corrélation de Pearson de la fréquence des indel parmi les différentes répétitions d’expériences. d La distribution de l’activité des gRNA pour les trois nucléases Cas9
Doench et al. ont développé un outil en ligne qui permet de concevoir des gRNA pour le knockout de gènes avec WT-SpCas925. Cet outil scanne une séquence codant pour un gène entier et classe tous les gRNA en fonction de leur activité et des effets hors cible. Nous avons utilisé cet outil pour concevoir des ARNg pour le criblage de la bibliothèque. Quatre ARNg les mieux classés, initiés par A ou G, ont été sélectionnés pour chaque gène (Fig. 2b). Nous avons également conçu des ARNg ciblant les microARN. Comme la longueur des séquences codantes des microARN est beaucoup plus courte que celle des régions codantes des gènes, nous avons généralement conçu trois ARNg pour chaque microARN. Un total de 80 263 oligonucléotides contenant des gRNA et les séquences cibles correspondantes (75 312 gRNA pour 19 037 gènes codants ; 4951 gRNA pour 1549 microARN) ont été synthétisés par microarray (Données supplémentaires 1). Les oligonucléotides ont été amplifiés par PCR et clonés dans les vecteurs lentivirus par assemblage de Gibson. L’analyse de la bibliothèque de plasmides par séquençage profond a révélé que le taux d’erreur (une lecture contenant une quelconque mutation était considérée comme une erreur) induit par la synthèse des oligonucléotides ou l’amplification par PCR au niveau de la région de la séquence cible de l’ARN guide était de 36,5 %. Cette bibliothèque de plasmides a été utilisée dans les expériences de criblage groupé suivantes pour profiler l’activité des ARNg pour WT-SpCas9, eSpCas9(1.1) et SpCas9-HF1.
La bibliothèque a été empaquetée dans des lentivirus et transduite dans des cellules HEK293T exprimant WT-SpCas9, eSpCas9(1.1) ou SpCas9-HF1 à une MOI de 0,3. Après 5 jours d’édition du génome, l’ADN génomique a été extrait, et les régions cibles intégrées ont été amplifiées par PCR pour un séquençage profond (Fig. 2b). Les mutations au niveau des régions de séquence cible de l’ARN guide peuvent être induites soit par les nucléases Cas9, soit par la construction de librairies. Si une mutation peut être trouvée dans la bibliothèque originale, elle a été considérée comme une mutation induite par la construction de la bibliothèque et exclue de l’analyse des indels. Les indels ont pu être détectés par séquençage profond sur les sites cibles intégrés (figure supplémentaire 6a). Nous avons obtenu des taux d’indels gRNA valides (nombre de lectures > 100) de 55 604 (couvrant 20 211 gènes), 58 167 (couvrant 20 315 gènes), et 56 888 (couvrant 20 270 gènes) pour WT-SpCas9, eSpCas9(1.1), et SpCas9-HF1, respectivement (Données supplémentaires 2). À notre connaissance, il s’agit des plus grands ensembles d’activité sur cible des gRNA rapportés jusqu’à présent dans les cellules de mammifères.
Le test de criblage a été répété expérimentalement deux fois, et deux répliques indépendantes ont montré un niveau élevé de corrélation pour le taux d’indel (R = 0,92 pour WT-SpCas9 ; R = 0,89 pour eSpCas9(1.1) ; R = 0,91 pour SpCas9-HF1, figure 2c). Le taux d’indel des gRNA individuels présente également une forte corrélation entre les trois nucléases Cas9 (Fig. 2c), ce qui indique que certaines caractéristiques de la séquence sont favorisées pour ces trois nucléases Cas9. La distribution des activités des gRNA varie remarquablement, de l’absence d’activité à un taux d’indel de 100 % pour ces trois nucléases Cas9 (Fig. 2d). WT-SpCas9 a montré une plus grande efficacité d’édition que eSpCas9(1.1) et SpCas9-HF1 dans notre criblage (Fig. 6b supplémentaire). Étant donné qu’il s’agissait de clones dérivés de cellules uniques pouvant influencer les efficacités, nous avons sélectionné cinq gRNA et les avons transfectés avec la nucléase Cas9 individuelle dans des cellules. WT-SpCas9 et SpCas9-HF1 ont montré une activité similaire, mais eSpCas9(1.1) a montré une activité plus faible au jour 3 (figure supplémentaire 6c), ce qui est cohérent avec les études précédentes15,16.
Il a été signalé que l’ADN plasmidique résiduel des procédures d’emballage viral peut contaminer les cellules transduites35, ce qui entraîne des inexactitudes potentielles dans la mesure de l’activité des gRNA. Nous avons conçu une paire d’amorces spécifiques du squelette des plasmides pour détecter les plasmides résiduels, et une paire d’amorces spécifiques de l’ADN génomique du lentivirus pour détecter à la fois les plasmides résiduels et les lentivirus intégrés dans le génome (figure supplémentaire 7a). Ces deux paires d’amorces ont montré une efficacité d’amplification similaire lorsque l’ADN plasmidique était utilisé comme modèle (figure supplémentaire 7b). L’ADN plasmidique résiduel a pu être détecté à la fois dans les virus non concentrés et dans les virus concentrés pendant le conditionnement du virus (figure supplémentaire 7b). Nous avons transduit des cellules HEK293T et extrait l’ADN génomique au jour 1 et au jour 5 après la transduction, respectivement. L’ADN plasmidique résiduel a pu être détecté dans les deux échantillons, mais les bandes PCR étaient très faibles au jour 5 avec l’amorce spécifique du squelette, ce qui indique que le plasmide résiduel s’est dégradé avec le temps (figure supplémentaire 7b). En revanche, des bandes très fortes ont pu être détectées avec les amorces spécifiques aux virus au jour 5, ce qui indique que de plus en plus de lentivirus se sont intégrés au génome. Ces résultats suggèrent que l’ADN plasmidique résiduel n’a eu qu’une influence minime sur le criblage de la bibliothèque.
Les caractéristiques de séquence associées à l’activité des gRNA
La caractérisation des caractéristiques de séquence associées à l’activité des gRNA est cruciale pour le développement d’outils de conception de gRNA. L’ensemble de données à grande échelle généré ici nous permet de mieux évaluer quelles caractéristiques ont le plus contribué à l’activité des gRNA. Des algorithmes tels que les arbres de régression boostés par gradient et la régression lasso ont été utilisés pour évaluer l’importance des caractéristiques36. Cependant, les arbres de régression boostés par gradient fournissent des scores d’importance de Gini qui ne reflètent que la valeur absolue de la contribution des caractéristiques, ce qui entraîne la perte d’informations concernant la direction de l’effet ; la régression lasso présente une faible capacité descriptive. Heureusement, un algorithme récemment développé SHAP (SHapley Additive exPlanation), une approche unifiée pour expliquer la sortie de tout modèle d’apprentissage automatique, peut potentiellement répondre à ces limitations37.
Nous avons connecté XGBoost avec SHAP (appelé Tree SHAP) pour évaluer l’importance de 1031 caractéristiques, y compris les caractéristiques identifiées par Doench et Wong et al.25,38, et toutes les accessibilités de position des caractéristiques de la structure secondaire des gRNA (Données supplémentaires 3, 4). Dans l’ensemble, les scores prédits ont été fortement influencés par la composition nucléotidique dépendant de la position pour trois nucléases Cas9 (Fig. 3a-c ; Données supplémentaires 4). Le nucléotide le plus favorisé était G en position 20 (G_20), le nucléotide immédiatement adjacent à la séquence PAM. D’autres caractéristiques importantes se chevauchent dans le top 20 pour trois nucléases : la température de fusion (Tm, favorisée), le nombre de dimères TT (défavorisé), C_18 (favorisé), l’énergie libre d’auto-pliement (favorisée) et G_14 (défavorisé).
Analyse de l’importance des caractéristiques associées à l’activité des gRNA par Tree SHAP. a-c 20% des caractéristiques les plus importantes identifiées par Tree SHAP pour WT-SpCas9, eSpCas9(1.1) et SpCas9-HF1, respectivement. Les nucléotides, ainsi que leur position, sont indiqués à gauche. GG_19 signifie le début du dimère GG à la position 19. Tm signifie température de fusion
Nous avons ensuite évalué la composition nucléotidique dépendante de la position des 25 % d’ARNg actifs les plus élevés par rapport aux 25 % d’ARNg actifs les plus bas. Les résultats ont révélé que G était généralement favorisé, et que T était généralement défavorisé (Fig. 4a-c ; Données supplémentaires 5). G_20 était fortement favorisé pour les trois nucléases Cas9, ce qui correspond à l’analyse SHAP de l’arbre. Les différences de préférence nucléotidique entre les variantes WT-SpCas9 et SpCas9 ont également été observées. Les différences de nucléotide favorisé en position 3 (C/G vs G), 9 (G vs C/G), 10 (A/G vs A/C), 14 (C vs A/T), 16 (A/C vs C), 17 (G vs A) et 18 (C/G vs C) ont été observées entre WT-SpCas9 et eSpCas9(1.1). Les différences de nucléotide favorisé en position 3 (C/G vs G), 5 (G vs C/T), 7 (C/G vs G), 9 (G vs C/G), 10 (A/G vs C), 12 (A/G vs A), 14 (C vs A/T), 17 (G vs A/C), et 18 (C/G vs C) ont été observées entre WT-SpCas9 et SpCas9-HF1. En outre, nous avons analysé la composition nucléotidique indépendante de la position avec les 20 % d’ARNg les plus actifs (Fig. 4d-f). Le contenu G (favorisé) et T (défavorisé) a fortement influencé l’activité des gRNA, tandis que le contenu A et C a légèrement influencé l’activité des gRNA pour trois nucléases Cas9, ce qui est cohérent avec les valeurs SHAP de l’arbre.
Influence de la composition nucléotidique sur l’activité des ARNg. a-c La composition nucléotidique dépendant de la position des 25 % d’ARNg actifs les plus élevés par rapport aux 25 % d’ARNg actifs les plus bas. Les barres montrent les scores log-odds de la fréquence des nucléotides pour chaque position. Les chiffres en dessous indiquent la position des nucléotides sur l’ADN cible. d-f L’association de chaque numéro de nucléotide avec l’activité du gRNA. La taille des cercles indiquait la fréquence des indel
Performance des algorithmes conventionnels
En plus de générer des ensembles de données sur l’activité des gRNA, un autre objectif de ce travail était de développer des outils de prédiction pour la conception de gRNA. Nous avons évalué les performances de quatre algorithmes conventionnels de prédiction de l’activité des gRNA, notamment la régression linéaire, la régression linéaire régularisée L2 (régression Ridge), la régression XGBoost et les modèles de perceptron multicouche (MLP) avec les ensembles de données générés dans cette étude. Pour éviter un ajustement excessif, nous avons séparé aléatoirement l’ensemble de données en deux sous-groupes, 85 % des données étant utilisées comme ensemble de données de formation pour former les modèles, et les 15 % restants étant utilisés pour tester la capacité de généralisation des modèles formés (Fig. 5a). Pour obtenir des performances optimales, les caractéristiques présentant des valeurs élevées de Tree SHAP ont été modélisées dans les algorithmes (données supplémentaires 4).
Performance des différents algorithmes pour la prédiction de l’activité des gRNA. a Schéma de l’ensemble des données et des algorithmes conventionnels. Quatre algorithmes conventionnels comprenant la régression linéaire, la régression ridge, la régression XGBoost et le MLP ont été construits, respectivement. Au total, 85 % de l’ensemble des données pertinentes ont été utilisées comme ensemble d’apprentissage, et les 15 % restants de l’ensemble des données dans chaque ensemble comme ensemble de test pour mesurer la capacité de généralisation de chaque modèle pour prédire des données non vues. b Schéma de l’ensemble des données et des algorithmes d’apprentissage profond. Des algorithmes d’apprentissage profond comprenant CNN, RNN et RNN + biofeature ont été construits, respectivement. Au total, 76,5 % de l’ensemble de données pertinentes ont été utilisées comme ensemble d’apprentissage, 8,5 % de l’ensemble de données de chaque ensemble ont été utilisées comme ensemble de validation et les 15 % restants de l’ensemble de données de chaque ensemble ont été utilisées comme ensemble de test pour mesurer la capacité de généralisation de chaque modèle pour prédire les données non vues. c-e Performance des différents algorithmes pour la prédiction de l’activité des ARNg révélée par la corrélation de Spearman pour WT-SpCas9, eSpCas9(1.1) et SpCas9-HF1, respectivement. Le graphique en barres montre la moyenne ± l’écart-type pour le coefficient de corrélation de Spearman entre les scores d’activité gRNA prédits et mesurés (n = 10)
Parmi les quatre algorithmes testés ici, MLP est le plus prédictif, avec des coefficients de corrélation de Spearman de 0.8416, 0,8457, et 0,8440 pour WT-SpCas9, eSpCas9(1.1), et SpCas9-HF1, respectivement (Fig. 5c-e ; Données supplémentaires 6-8). XGBoost est le deuxième meilleur prédicteur, avec des coefficients de corrélation de Spearman de 0,8454, 0,8310 et 0,8184 pour WT-SpCas9, eSpCas9(1.1) et SpCas9-HF1, respectivement. La régression linéaire et la régression ridge ont également donné de bons résultats, mais avec un score de corrélation relativement plus faible. Nous avons également essayé la régression lasso et la régression SVM. Mais le coefficient de pénalité de la régression lasso est presque nul, ce qui la rend équivalente au modèle linéaire. La régression SVM n’a pas réussi à terminer le benchmarking sur l’échelle actuelle de l’ensemble de données dans un délai raisonnable (3 semaines). Ils ont été abandonnés dans la comparaison finale.
Performance des algorithmes d’apprentissage profond
Des études récentes ont montré que deux algorithmes basés sur l’apprentissage profond, le réseau neuronal convolutif (CNN) et le réseau neuronal récurrent (RNN), sont des outils puissants pour l’analyse liée aux séquences d’ADN/protéines39,40,41,42,43. Ils permettent d’obtenir automatiquement des caractéristiques utiles à partir de séquences brutes d’ADN/protéines, sans qu’il soit nécessaire de recourir à l’ingénierie des caractéristiques. CNN a été utilisé pour prédire l’activité gRNA pour Cpf1 et WT-SpCas926,27, tandis que RNN n’a pas été utilisé pour la prédiction de l’activité gRNA jusqu’à présent. Nous avons entraîné à la fois CNN et RNN pour la prédiction de l’activité gRNA. Pour éviter l’ajustement excessif, nous avons séparé aléatoirement l’ensemble de données en trois sous-groupes avec 76,5 % des données utilisées comme ensemble de données de formation pour former les modèles, 8,5 % utilisés comme ensemble de données de validation, et les 15 % restants retenus utilisés pour tester la capacité de généralisation des modèles formés (Fig.
RNN a surpassé CNN et les autres algorithmes pour la prédiction de l’activité gRNA avec des coefficients de corrélation de Spearman de 0,8555, 0,8491 et 0,8512 pour WT-SpCas9, eSpCas9(1.1) et SpCas9-HF1, respectivement (Fig. 5c-e ; Données supplémentaires 6-8). Le CNN a obtenu des performances similaires à celles de XGBoost avec des coefficients de corrélation de Spearman de 0,8455, 0,8313 et 0,8343 pour WT-SpCas9, eSpCas9(1.1) et SpCas9-HF1, respectivement.
Un modèle intégré améliore le pouvoir prédictif
L’épine dorsale des algorithmes d’apprentissage profond (CNN ou RNN) ne peut exploiter que la composition des k-mer ou ses dépendances39,44. Des études récentes sur la prédiction liée aux protéines ont montré que la capacité de prédiction des modèles d’apprentissage profond pourrait être stimulée par l’ajout d’autres caractéristiques, telles que le poids moléculaire, l’hydrophobie et la charge absolue, qui ne pourraient pas être obtenues automatiquement par les modèles d’apprentissage profond45,46. Dans notre travail, les caractéristiques indirectes de la séquence, y compris les accessibilités de position de la structure secondaire, la tige-boucle de la structure secondaire, la température de fusion et la teneur en GC sont fortement associées à l’activité de l’ARNg (données supplémentaires 4), mais elles n’ont pas pu être obtenues par apprentissage profond. Considérant que le modèle RNN a obtenu la meilleure performance de tous les algorithmes, nous avons donc combiné ces caractéristiques biologiques avec RNN pour la prédiction de l’activité gRNA. Il est intéressant de noter que l’ajout de ces caractéristiques au RNN a augmenté la puissance de prédiction, avec des coefficients de corrélation de Spearman de 0,8670, 0,8624 et 0,8603 pour WT-SpCas9, eSpCas9(1.1) et SpCas9-HF1, respectivement (Fig. 5c-e). Par conséquent, le RNN intégré aux caractéristiques biologiques (ci-après dénommé RNN + Biofeature) a été utilisé comme modèle final pour la prédiction de l’activité des gRNA (Fig. 6).
Modèle DeepHF pour la prédiction de l’activité des gRNA. La séquence originale du gRNA est d’abord encodée, puis intégrée pour obtenir une nouvelle représentation. Cette nouvelle représentation est ensuite traitée par un BiLSTM pour obtenir la représentation finale qui est ensuite concaténée avec des caractéristiques fabriquées à la main pour servir d’entrée à la couche entièrement connectée pour être transformée de manière non linéaire. Enfin, une transformation linéaire est effectuée pour obtenir le score de prédiction
Pour tester davantage les performances des sept modèles utilisés dans cette étude, nous avons généré une liste de taux d’indel gRNA pour les sites endogènes (85 sites pour WT-SpCas9, 81 sites pour eSpCas9(1.1) et 82 sites pour SpCas9-HF1) (Données supplémentaires 9). Les sept modèles fonctionnaient très bien, mais la régression linéaire et la régression ridge étaient moins prédictives d’après les mesures de corrélation de Spearman (figures supplémentaires 8-10). En raison de l’ensemble de données limité ici, nous n’avons pas pu conclure quel algorithme était statistiquement meilleur que les autres algorithmes sur les sites endogènes. Nous avons également étudié la corrélation de la fréquence des indel sur les cibles synthétiques avec celle des cibles endogènes correspondantes. La corrélation de Spearman est de 0,722, 0,767 et 0,730 pour WT-SpCas9, eSpCas9(1.1) et SpCas9-HF1, respectivement (figure supplémentaire 11a-c).
Comparaison du modèle RNN + Biofeature avec les modèles existants
Il existe plusieurs ensembles de données publiquement disponibles sur l’efficacité des gRNA pour WT-SpCas9, ce qui nous permet de comparer les performances du modèle RNN + Biofeature pour WT-SpCas9 (ci-après dénommé DeepWt) avec les modèles de prédiction existants. Nous avons testé DeepWt par rapport à 18 ensembles de données endogènes collectées par Haeussler et al.47, et les coefficients de corrélation de Spearman obtenus varient de 0,129 à 0,594 (Données supplémentaires 10). Il a été signalé que le modèle de prédiction dépend fortement du fait que l’ARNg soit exprimé à partir d’un promoteur U6 dans les cellules ou d’un promoteur T7 in vitro47,48. Par conséquent, la stratégie d’apprentissage par transfert a été utilisée pour améliorer la capacité de prédiction de notre modèle dans différentes conditions d’expression. Pour l’expression du promoteur U6, les performances ont été améliorées par un réglage fin de la dernière couche cachée de DeepWt avec l’ensemble de données XuKBM (données supplémentaires 10). Notre modèle final optimisé, appelé DeepWt_U6, a surpassé les performances des sept autres outils de conception de gRNA populaires (Fig. 7a). Notamment, puisque notre bibliothèque de gRNA a été conçue, en partie, par RuleSet2 qui a été développé par les ensembles de données Doench, cette comparaison peut avoir un biais pour DeepWt_U6 sur les ensembles de données Doench. Pour l’expression du promoteur T7, nous avons développé un autre modèle appelé DeepWt_T7 en affinant l’algorithme RNN + biofeature avec l’ensemble de données Moreno-Mateos2015. Ce modèle a surpassé les autres modèles pour la conception de gRNAs exprimés in vitro (Fig. 7b).
Carte thermique des coefficients de corrélation de rang de Spearman entre les scores d’efficacité et les ensembles de données. a ARNg transcrits dans les cellules à partir d’un promoteur U6. b ARNg transcrits in vitro à partir d’un promoteur T7. Le type de cellule, le nombre de gRNA et l’espèce sont indiqués à gauche. Les scores sont indiqués sur l’axe horizontal, les ensembles de données sur l’axe vertical. Les corrélations d’un algorithme par rapport à son propre ensemble de données d’entraînement sont indiquées en gris car elles sont susceptibles d’être surestimées en raison d’un ajustement excessif. Les scores les plus élevés sont en gras
Il a été signalé que l’intégration des métriques d’accessibilité du site cible dans le modèle pourrait améliorer le pouvoir de prédiction26,27. Pour WT-SpCas9, nous avons affiné le modèle DeepWt_U6 avec des données de DNase I de la lignée cellulaire KBM-7 obtenues à partir d’ENCODE, ce qui a donné le modèle DeepWt_Chromatin. Les ensembles de données Wang/Xu HL60 et Hart Hct116-2 Lib 1 ainsi que les données DNase I correspondantes ont été utilisés pour tester les performances de DeepWt_Chromatin. Cependant, les scores de corrélation de Spearman n’ont pas été améliorés (Données supplémentaires 11). Nous avons également récupéré les données de DNase I des cellules HEK293T de la base de données ENCODE et avons testé si l’intégration des métriques dans les modèles pouvait améliorer le pouvoir de prédiction pour eSpCas9(1.1) et SpCas9-HF1. Les données relatives à la DNase I ont été traitées selon la méthode décrite par Kim et al.26. Cependant, l’incorporation de ces données n’a pas pu augmenter significativement la capacité de prédiction (validation shuffled décuplée) pour eSpCas9(1.1) et SpCas9-HF1 (Données supplémentaires 12).
Les contributions nucléotidiques révélées par Deep SHAP
En plus de la précision de prédiction, nous sommes également intéressés par la compréhension des mécanismes du modèle d’apprentissage profond. Lundberg et Lee49 ont développé un algorithme appelé Deep SHAP, qui est un algorithme d’approximation rapide des valeurs SHAP dans les modèles d’apprentissage profond. Nous avons utilisé Deep SHAP pour estimer la contribution des nucléotides dépendant de la position dans le modèle d’apprentissage profond. La contribution de chaque nucléotide dépendant de la position à l’activité du gRNA a été calculée à partir de la valeur moyenne de cette position dans tous les gRNA d’entraînement. Pour rendre la contribution des nucléotides comparable entre WT-SpCas9, eSpCas9(1.1) et SpCas9-HF1, les valeurs SHAP profondes ont été rééchelonnées par le Z-score (c’est-à-dire la normalisation). Le nucléotide avec un Z-score supérieur à 1 ou inférieur à -1 a été considéré comme ayant une contribution significative à l’activité du gRNA (figure supplémentaire 12a).
Il y avait 16, 21 et 26 nucléotides significatifs dans WT-SpCas9, eSpCas9(1.1) et SpCas9-HF1, respectivement. Le résultat a révélé que G avait typiquement une contribution positive et T typiquement une contribution négative, en accord avec l’observation précédente que Cas9 se lie préférentiellement aux ARNg contenant des purines mais pas des pyrimidines50. De plus, des Ts multiples dans l’espaceur provoquaient une faible expression de l’ARNg51. Nous avons constaté que la plupart des nucléotides significatifs avaient la même direction de contribution à l’activité gRNA dans les trois nucléases Cas9. Conformément à plusieurs rapports précédents38,52, les nucléotides de la position 20 étaient les plus influents, G_20 ayant une forte contribution positive et C_20/T_20 une forte contribution négative. Par rapport à WT-SpCas9, les motifs spécifiques à eSpCas9(1.1) comprennent A_15 (favorisé), A_17 (favorisé), G_6-8 (favorisé), G_14 (défavorisé), G_16 (défavorisé), T_6 (défavorisé), T_11-12 (défavorisé) ; Les motifs spécifiques à SpCas9-HF1 comprenaient A_11 (favorisé), A_13 (défavorisé), A_14-15 (favorisé), A_19-20 (défavorisé), G_6-7 (favorisé), G_14 (défavorisé), G_16 (défavorisé), T_4 (défavorisé), T_6 (défavorisé) et T_11 (défavorisé) (figure supplémentaire. 12a).
En outre, la différence de Z-score entre les variantes de Cas9 et WT-SpCas9 a été calculée pour évaluer les changements dans les contributions nucléotidiques. La différence supérieure à 1 ou inférieure à -1 a été considérée comme ayant un changement significatif (figure supplémentaire 12b). Plusieurs différences dans le sens de la contribution des nucléotides significatifs ont été observées. Plus précisément, A_13 a contribué négativement à SpCas9-HF1 mais pas à WT-SpCas9, G_17 a contribué positivement à WT-SpCas9 mais pas à eSpCas9(1.1) ni à HF1-SpCas9, G_18 a contribué positivement à WT-SpCas9 mais négativement à eSpCas9(1.1) (également à HF1-SpCas9, mais de manière non significative).
La contribution des nucléotides répétitifs aux activités des gRNA a également été étudiée sur la base de la somme des valeurs SHAP profondes. Pour WT-SpCas9, une étude précédente a montré qu’un tronçon de nucléotides identiques adjacents (nucléotides répétitifs, tels que GGGG ou TTTT) pourrait être associé à une mauvaise efficacité d’un sgRNA38. Cependant, ils n’ont pas pris en compte les effets de position. Par conséquent, nous avons calculé la moyenne des valeurs SHAP profondes des quatre nucléotides répétitifs comprenant AAAA, CCCC, GGGG, et TTTT de la position 1 à 18 (Données supplémentaires 13). Nos résultats ont démontré que les nucléotides répétitifs diminuent généralement l’efficacité des indel, ce qui est cohérent avec l’étude de Wong et al38. Cependant, la contribution positive des nucléotides répétitifs à certaines positions a été observée, notamment GGGG à partir de la position 1-5 et 16-18, AAAA à partir de la position 14, CCCC à partir de la position 1-2 et 15-18 (Données supplémentaires 13). La contribution positive des nucléotides répétés à certaines positions a également été observée pour eSpCas9 (1.1) et SpCas9-HF1. Par exemple, AAAA à partir de la position 14-15 a contribué positivement aux activités gRNA pour les deux nucléases Cas9 (Supplementary Data 13).
La corrélation entre la fréquence d’indel et le phénotype
Dans cette étude, nous avons utilisé le taux d’indel comme label d’activité gRNA, ce qui n’est pas égal à l’efficacité réelle du knockout génique. Nous avons testé la corrélation entre la fréquence des indels et la perturbation réelle des gènes à l’aide d’un test basé sur les protéines. Nous avons conçu un total de neuf ARNg ciblant SIRT1, SIRT2 et SIRT6, avec trois ARNg pour chacun. Les gRNA et les nucléases Cas9 ont été introduits dans des cellules HEK293T à l’aide d’un vecteur épisomal qui permet l’édition du génome à long terme5. La fréquence des indels et l’expression des protéines ont été analysées au jour 9 après la transfection. Les résultats ont révélé que les fréquences d’indel présentaient une bonne corrélation avec l’expression des protéines (r = 0,82 ; figures 13a et b supplémentaires). En outre, la corrélation entre la fréquence des indels et la perturbation effective du gène a été testée par un test de rapporteur à la luciférase. Nous avons conçu un total de 11 ARNg ciblant le gène de la luciférase. Cinq jours après l’édition du génome, la fréquence des indels et l’activité de la luciférase ont été analysées. Les résultats ont révélé que les fréquences d’indel avaient une bonne corrélation avec l’activité de la luciférase (r = 0,70, figure supplémentaire 13c-e).
La corrélation entre l’efficacité on-target et off-target
Le but de ce travail est de concevoir des gRNA avec une meilleure activité on-target, mais ces gRNA peuvent avoir une plus grande tolérance aux mésappariements et donc induire des mutations off-target plus élevées. Pour tester cette hypothèse, nous avons sélectionné trois sgRNA ayant une activité différente et conçu des paires guide-cible avec des doubles mésappariements (figure supplémentaire 14a). Le clivage hors cible s’est produit efficacement avec les séquences cibles contenant des mésappariements en position 1-8 pour WT-SpCas9 (figure supplémentaire 14b). Comme prévu, une activité on-target plus élevée entraîne généralement une activité off-target plus élevée. Pour eSpCas9(1.1) et SpCas9-HF1, cependant, le clivage hors cible était au niveau de fond, sauf pour les mésappariements en position 7-8, où le site1-gRNA a montré un clivage hors cible élevé (figure supplémentaire 14b). Une tolérance similaire aux mésappariements a également été observée par Slaymaker et al.15 .. Le site2-gRNA avait une activité comparable à celle du site1-gRNA, mais son clivage hors cible était au niveau de fond, ce qui indique que la tolérance au mismatch dépend des séquences gRNA.
Service en ligne
Nous avons finalement développé un outil en ligne appelé DeepHF (Deep learning for High-Fidelity Cas9) basé sur le modèle RNN + biofeature pour la conception de gRNA pour WT-SpCas9, eSpCas9(1.1) et SpCas9-HF1. L’outil en ligne contient trois modules fonctionnels, à savoir le module de prédiction, le module des gRNA vérifiés et le module de conception. Le module de prédiction permet aux utilisateurs d’obtenir les activités prédites pour tous les gRNA avec une séquence d’ADN d’entrée. Le module des gRNA vérifiés fournit tous les taux d’indel des gRNA générés dans cette étude. Le module de conception fournit les ARNg qui conviennent à l’élimination des gènes avec eSpCas9(1.1) et SpCas9-HF1 dans les cellules humaines. Dans ce module, les gRNA ont été choisis parmi les transcrits communs de chaque gène (Genome Reference Consortium Human Build 38). Les informations hors cible (1-3 mésappariements considérés comme hors cible) et l’emplacement du ciblage (dans les 5-65% de la séquence codante) ont également été annotés. Les utilisateurs peuvent obtenir les ARNg prédéfinis en saisissant un ID de gène ou un symbole de gène. Le site web est disponible gratuitement à l’adresse http://www.DeepHF.com/.