Un guide étape par étape de l’analyse des données de séquençage de l’ADN

Le Dr Javier Quilez Oliete, consultant expérimenté en bioinformatique sur Kolabtree, fournit un guide complet de l’analyse des données de séquençage de l’ADN, y compris les outils et les logiciels utilisés pour lire les données.

Introduction

L’acide désoxyribonucléique (ADN) est la molécule qui porte la plupart des informations génétiques d’un organisme. (Dans certains types de virus, l’information génétique est portée par l’acide ribonucléique (ARN).) Les nucléotides (conventionnellement représentés par les lettres A, C, G ou T) sont les unités de base des molécules d’ADN. Conceptuellement, le séquençage de l’ADN est le processus de lecture des nucléotides qui composent une molécule d’ADN (par exemple, « GCAAACCAAT » est une chaîne d’ADN de 10 nucléotides). Les technologies de séquençage actuelles produisent des millions de lectures d’ADN en un temps raisonnable et à un coût relativement faible. À titre de référence, le coût du séquençage d’un génome humain – un génome est l’ensemble des molécules d’ADN d’un organisme – a franchi la barre des 100 dollars et peut être réalisé en quelques jours. Cela contraste avec la première initiative de séquençage du génome humain, qui a été achevée en une décennie et a coûté environ 2,7 milliards de dollars.

Cette capacité de séquencer l’ADN à haut débit et à faible coût a permis le développement d’un nombre croissant de méthodes et d’applications basées sur le séquençage. Par exemple, le séquençage de génomes entiers ou de leurs régions codant pour les protéines (deux approches connues respectivement sous le nom de séquençage du génome entier et de l’exome) chez des individus malades et sains peut faire allusion à des altérations de l’ADN causant des maladies. De même, le séquençage de l’ARN qui est transcrit à partir de l’ADN – une technique connue sous le nom de séquençage de l’ARN – est utilisé pour quantifier l’activité des gènes et la façon dont elle change dans différentes conditions (par exemple, sans traitement ou avec traitement). D’autre part, les méthodes de séquençage par capture de la conformation des chromosomes détectent les interactions entre les molécules d’ADN proches et aident ainsi à déterminer la distribution spatiale des chromosomes dans la cellule.

Le point commun entre ces applications et d’autres applications du séquençage de l’ADN est la génération de jeux de données de l’ordre du gigaoctet et comprenant des millions de séquences lues. Par conséquent, pour donner un sens aux expériences de séquençage à haut débit (HTS), il faut disposer d’importantes capacités d’analyse des données. Heureusement, des outils informatiques et statistiques spécialisés et des flux d’analyse relativement standard existent pour la plupart des types de données HTS. Si certaines des étapes d’analyse (initiales) sont communes à la plupart des types de données de séquençage, l’analyse en aval dépendra du type de données et/ou de l’objectif final de l’analyse. Je fournis ci-dessous une introduction aux étapes fondamentales de l’analyse des données HTS et je fais référence à des outils populaires.

Certaines des sections ci-dessous sont axées sur l’analyse des données générées par les technologies de séquençage à lecture courte (principalement Illumina), car celles-ci ont historiquement dominé le marché HTS. Cependant, des technologies plus récentes qui génèrent des lectures plus longues (par exemple Oxford Nanopore Technologies, PacBio) gagnent rapidement du terrain. Comme le séquençage à longue lecture présente certaines particularités (par exemple, des taux d’erreur plus élevés), des outils spécifiques sont développés pour l’analyse de ce type de données.

Contrôle de qualité (QC) des lectures brutes

L’analyste avide commencera l’analyse à partir de fichiers FASTQ ; le format FASTQ est depuis longtemps la norme pour stocker les données de séquençage à lecture courte. En substance, les fichiers FASTQ contiennent la séquence de nucléotides et la qualité d’appel par base pour des millions de lectures. Bien que la taille du fichier dépende du nombre réel de lectures, les fichiers FASTQ sont généralement volumineux (de l’ordre des mégaoctets et des gigaoctets) et compressés. Il est à noter que la plupart des outils qui utilisent les fichiers FASTQ en entrée peuvent les traiter en format compressé. Ainsi, afin d’économiser de l’espace disque, il est recommandé de ne pas les décompresser. Par convention, j’assimilerai ici un fichier FASTQ à un échantillon de séquençage.

FastQC est probablement l’outil le plus populaire pour effectuer le CQ des lectures brutes. Il peut être exécuté via une interface visuelle ou par programme. Si la première option peut être plus pratique pour les utilisateurs qui ne se sentent pas à l’aise avec l’environnement en ligne de commande, la seconde offre une évolutivité et une reproductibilité incomparables (pensez à la façon dont il peut être fastidieux et sujet aux erreurs d’exécuter manuellement l’outil pour des dizaines de fichiers). Quoi qu’il en soit, le résultat principal de FastQC est un fichier HTML présentant des statistiques sommaires clés sur la qualité globale des lectures de séquençage brutes d’un échantillon donné. Inspecter des dizaines de rapports FastQC un par un est fastidieux et complique la comparaison entre les échantillons. Par conséquent, vous pourriez vouloir utiliser MultiQC, qui regroupe les rapports HTML de FastQC (ainsi que d’autres outils utilisés en aval, par exemple le trimming d’adaptateur, l’alignement) dans un seul rapport.

Les informations de QC sont destinées à permettre à l’utilisateur de juger si les échantillons ont une bonne qualité et peuvent donc être utilisés pour les étapes suivantes ou s’ils doivent être écartés. Malheureusement, il n’existe pas de seuil consensuel basé sur les métriques FastQC pour classer les échantillons comme étant de bonne ou de mauvaise qualité. L’approche que j’utilise est la suivante. Je m’attends à ce que tous les échantillons qui sont passés par la même procédure (par exemple, l’extraction d’ADN, la préparation de la bibliothèque) aient des statistiques de qualité similaires et une majorité de drapeaux « pass ». Si certains échantillons ont une qualité inférieure à la moyenne, je les utiliserai quand même dans l’analyse en aval en gardant cela à l’esprit. D’un autre côté, si tous les échantillons de l’expérience obtiennent systématiquement des drapeaux « warning » ou « fail » dans plusieurs métriques (voir cet exemple), je soupçonne que quelque chose n’a pas fonctionné dans l’expérience (par exemple, une mauvaise qualité d’ADN, la préparation de la bibliothèque, etc.) et je recommande de la répéter.

Taillage des lectures

La QC des lectures brutes aide à identifier les échantillons problématiques, mais elle n’améliore pas la qualité réelle des lectures. Pour ce faire, nous devons découper les lectures pour éliminer les séquences techniques et les extrémités de mauvaise qualité.

Les séquences techniques sont des restes de la procédure expérimentale (par exemple, les adaptateurs de séquençage). Si de telles séquences sont adjacentes à la véritable séquence de la lecture, l’alignement (voir ci-dessous) peut cartographier les lectures à la mauvaise position dans le génome ou diminuer la confiance dans un alignement donné. Outre les séquences techniques, on peut également vouloir éliminer les séquences d’origine biologique si elles sont très présentes parmi les lectures. Par exemple, des procédures de préparation de l’ADN non optimales peuvent laisser une forte proportion d’ARN ribosomal (ARNr) converti en ADN dans l’échantillon. À moins que ce type d’acide nucléique ne soit la cible de l’expérience de séquençage, le fait de conserver les lectures dérivées de l’ARNr ne fera qu’augmenter la charge de calcul des étapes en aval et risque de fausser les résultats. À noter que si les niveaux de séquences techniques, d’ARNr ou d’autres contaminants sont très élevés, ce qui aura probablement déjà été mis en évidence par le CQ, vous pouvez vouloir rejeter l’ensemble de l’échantillon de séquençage.

Dans le séquençage à lecture courte, la séquence d’ADN est déterminée un nucléotide à la fois (techniquement, un nucléotide à chaque cycle de séquençage). En d’autres termes, le nombre de cycles de séquençage détermine la longueur de lecture. Un problème connu des méthodes de séquençage HTS est la diminution de la précision avec laquelle les nucléotides sont déterminés à mesure que les cycles de séquençage s’accumulent. Cela se traduit par une diminution globale de la qualité d’appel par base, en particulier vers la fin de la lecture. Comme cela se produit avec les séquences techniques, essayer d’aligner des lectures qui contiennent des extrémités de faible qualité peut conduire à un mauvais placement ou à une mauvaise qualité de cartographie.

Pour éliminer les séquences techniques/contaminantes et les extrémités de faible qualité, des outils de rognage de lecture comme Trimmomatic et Cutadapt existent et sont largement utilisés. Essentiellement, ces outils supprimeront les séquences techniques (disponibles en interne et/ou fournies par l’utilisateur) et rogneront les lectures en fonction de leur qualité tout en maximisant la longueur des lectures. Les lectures qui restent trop courtes après le découpage sont écartées (les lectures excessivement courtes, par exemple <36 nucléotides, compliquent l’étape d’alignement car elles sont susceptibles de correspondre à plusieurs sites dans le génome). Vous pouvez vouloir regarder le pourcentage de lectures qui survivent à l’élagage, car un taux élevé de lectures rejetées est probablement un signe de données de mauvaise qualité.

Enfin, je réexécute généralement FastQC sur les lectures rognées pour vérifier que cette étape a été efficace et a systématiquement amélioré les métriques de QC.

Alignement

A quelques exceptions près (par exemple l’assemblage de novo), l’alignement (également appelé cartographie) est généralement l’étape suivante pour la plupart des types de données et des applications HTS. L’alignement des lectures consiste à déterminer la position dans le génome d’où provient la séquence de la lecture (typiquement exprimée en tant que chromosome:start-end). Par conséquent, à cette étape, nous avons besoin de l’utilisation d’une séquence de référence pour aligner/mapper les lectures.

Le choix de la séquence de référence sera déterminé par de multiples facteurs. Tout d’abord, l’espèce dont provient l’ADN séquencé. Si le nombre d’espèces pour lesquelles une séquence de référence de haute qualité est disponible augmente, ce n’est peut-être toujours pas le cas pour certains organismes moins étudiés. Dans ces cas, vous pouvez aligner les lectures sur une espèce proche sur le plan évolutif pour laquelle un génome de référence est disponible. Par exemple, comme il n’existe pas de séquence de référence pour le génome du coyote, nous pouvons utiliser celle du chien, espèce étroitement apparentée, pour l’alignement des lectures. De même, nous pouvons toujours vouloir aligner nos lectures sur une espèce étroitement apparentée pour laquelle il existe une séquence de référence de meilleure qualité. Par exemple, si le génome du gibbon a été publié, celui-ci est fragmenté en milliers de fragments qui ne récapitulent pas entièrement l’organisation de ce génome en dizaines de chromosomes ; dans ce cas, effectuer l’alignement en utilisant la séquence de référence humaine peut être bénéfique.

Un autre facteur à prendre en compte est la version de l’assemblage de la séquence de référence, puisque de nouvelles versions sont publiées au fur et à mesure que la séquence est mise à jour et améliorée. Il est important de noter que les coordonnées d’un alignement donné peuvent varier selon les versions. Par exemple, de multiples versions du génome humain peuvent être trouvées dans le UCSC Genome Browser. Quelle que soit l’espèce, je favorise fortement la migration vers la version la plus récente de l’assemblage dès qu’elle est disponible. Cela peut causer quelques nuisances pendant la transition, car les résultats déjà existants seront relatifs aux anciennes versions, mais c’est payant à long terme.

A part cela, le type de données de séquençage a également de l’importance. Les lectures générées par les protocoles DNA-seq, ChIP-seq ou Hi-C seront alignées sur la séquence de référence du génome. D’autre part, comme l’ARN transcrit à partir de l’ADN est ensuite transformé en ARNm (c’est-à-dire que les introns sont supprimés), de nombreuses lectures d’ARN-seq ne pourront pas être alignées sur une séquence de référence du génome. Au lieu de cela, nous devons soit les aligner sur des séquences de référence du transcriptome, soit utiliser des aligneurs sensibles au fractionnement (voir ci-dessous) lorsque nous utilisons la séquence du génome comme référence. Le choix de la source pour l’annotation de la séquence de référence, c’est-à-dire la base de données contenant les coordonnées des gènes, des transcriptions, des centromères, etc. est lié à cette question. J’utilise généralement l’annotation GENCODE car elle combine une annotation complète des gènes et des séquences de transcription.

Une longue liste d’outils d’alignement de séquences à lecture courte a été développée (voir la section Alignement de séquences à lecture courte ici). Les passer en revue dépasse la portée de cet article (des détails sur les algorithmes derrière ces outils peuvent être trouvés ici). Selon mon expérience, les plus populaires sont Bowtie2, BWA, HISAT2, Minimap2, STAR et TopHat. Ma recommandation est que vous choisissiez votre aligneur en tenant compte de facteurs clés tels que le type de données HTS et l’application, ainsi que l’acceptation par la communauté, la qualité de la documentation et le nombre d’utilisateurs. Par exemple, on a besoin d’aligneurs comme STAR ou Bowtie2 qui sont conscients des jonctions exon-exon lors de la cartographie de l’ARN-seq au génome.

Le point commun à la plupart des mappeurs est la nécessité d’indexer la séquence utilisée comme référence avant que l’alignement réel ait lieu. Cette étape peut prendre du temps mais elle ne doit être effectuée qu’une seule fois pour chaque séquence de référence. La plupart des mappeurs stockent les alignements dans des fichiers SAM/BAM, qui suivent le format SAM/BAM (les fichiers BAM sont des versions binaires des fichiers SAM). L’alignement fait partie des étapes les plus longues et les plus complexes de l’analyse des données de séquençage et les fichiers SAM/BAM sont lourds (de l’ordre du gigaoctet). Il est donc important de s’assurer que vous disposez des ressources nécessaires (voir la dernière section ci-dessous) pour exécuter l’alignement dans un temps raisonnable et stocker les résultats. De même, en raison de la taille et du format binaire des fichiers BAM, évitez de les ouvrir avec des éditeurs de texte ; utilisez plutôt des commandes Unix ou des outils dédiés comme SAMtools.

D’après les alignements

Je dirais qu’il n’y a pas d’étape commune claire après l’alignement, car c’est à ce stade que chaque type de données HTS et chaque application peuvent différer.

Une analyse en aval courante pour les données ADN-seq est l’appel de variants, c’est-à-dire l’identification des positions dans le génome qui varient par rapport à la référence du génome et entre les individus. Un cadre d’analyse populaire pour cette application est GATK pour le polymorphisme nucléotidique simple (SNP) ou les petites insertions/délétions (indels) (figure 2). Les variants comprenant de plus gros morceaux d’ADN (également appelés variants structurels) nécessitent des méthodes d’appel dédiées (voir cet article pour une comparaison complète). Comme pour les aligneurs, je conseille de choisir le bon outil en tenant compte de facteurs clés comme le type de variants (SNP, indel ou variants structurels), l’acceptation par la communauté, la qualité de la documentation et le nombre d’utilisateurs.

L’application la plus fréquente de RNA-seq est probablement la quantification de l’expression génique. Historiquement, les lectures devaient être alignées sur la séquence de référence, puis le nombre de lectures alignées sur un gène ou un transcrit donné était utilisé comme un proxy pour quantifier ses niveaux d’expression. Cette approche alignement+quantification est réalisée par des outils tels que Cufflinks, RSEM ou featureCounts. Cependant, l’approche scuh a été de plus en plus dépassée par des méthodes plus récentes mises en œuvre dans des logiciels comme Kallisto et Salmon. Conceptuellement, avec ces outils, la séquence complète d’une lecture n’a pas besoin d’être alignée sur la séquence de référence. Au contraire, il suffit d’aligner suffisamment de nucléotides pour être sûr qu’une lecture provient d’un transcrit donné. En d’autres termes, l’approche alignement+quantification est réduite à une seule étape. Cette approche est connue sous le nom de pseudo-mapping et augmente considérablement la vitesse de quantification de l’expression génique. D’un autre côté, gardez à l’esprit que le pseudo-mapping ne conviendra pas aux applications pour lesquelles l’alignement complet est nécessaire (par exemple, l’appel de variants à partir de données RNA-seq).

Un autre exemple des différences dans les étapes d’analyse en aval et les outils requis à travers les applications basées sur le séquençage est ChIP-seq. Les lectures générées par cette technique seront utilisées pour l’appel de pic, qui consiste à détecter les régions du génome présentant un excès significatif de lectures qui indique où la protéine cible est liée. Plusieurs peak callers existent et cette publication les passe en revue. Comme dernier exemple, je mentionnerai les données Hi-C, dans lesquelles les alignements sont utilisés comme entrée pour des outils qui déterminent les matrices d’interaction et, à partir de celles-ci, les caractéristiques 3D du génome. Commenter tous les essais basés sur le séquençage dépasse le cadre de cet article (pour une liste relativement complète, voir cet article).

Avant de commencer…

La partie restante de cet article aborde des aspects qui peuvent ne pas être strictement considérés comme des étapes de l’analyse des données HTS et qui sont largement ignorés. En revanche, je soutiens qu’il est capital que vous réfléchissiez aux questions posées dans le tableau 1 avant de commencer à analyser les données HTS (ou tout type de données en fait), et j’ai écrit sur ces sujets ici et ici.

Tableau 1

Réfléchissez-y	Action proposée
Avez-vous toutes les informations de votre échantillon nécessaires à l’analyse ?	Collectez systématiquement les métadonnées des expériences
Pourrez-vous identifier sans équivoque votre échantillon ?	Etablissez un système pour attribuer à chaque échantillon un identifiant unique
Où se trouveront les données et les résultats ?	Organisation structurée et hiérarchique des données
Pourrez-vous traiter plusieurs échantillons de façon transparente ?	Échelle, parallélisation, configuration automatique et modularité du code
Pourrez-vous, ou quelqu’un d’autre, reproduire les résultats ?	Documentez votre code et vos procédures !

Comme mentionné ci-dessus, les données brutes HTS et certains des fichiers générés lors de leur analyse sont de l’ordre du gigaoctet, il n’est donc pas exceptionnel qu’un projet comprenant des dizaines d’échantillons nécessite des téraoctets de stockage. En outre, certaines étapes de l’analyse des données HTS sont gourmandes en ressources informatiques (par exemple, l’alignement). Cependant, l’infrastructure de stockage et de calcul requise pour l’analyse des données HTS est un élément important qui est souvent négligé ou non discuté. Par exemple, dans le cadre d’une analyse récente, nous avons examiné des dizaines d’articles publiés effectuant une analyse d’association à l’échelle du phénome (PheWAS). Les PheWAS modernes analysent 100 à 1 000 variantes génétiques et phénotypes, ce qui nécessite un stockage de données et une puissance de calcul importants. Et pourtant, pratiquement aucun des articles que nous avons examinés n’a commenté l’infrastructure nécessaire à l’analyse PheWAS. Il n’est pas surprenant que ma recommandation soit de planifier dès le départ les exigences de stockage et de calcul auxquelles vous serez confronté et de les partager avec la communauté.

Vous avez besoin d’aide pour analyser les données de séquençage de l’ADN ? Entrez en contact avec des bioinformaticiens indépendants et des experts en génomique sur Kolabtree.