TP module 8 bis

1 Objectifs

L’objectif du TP, à travers l’exemple de l’assemblage d’un souche bactérienne séquencée par HTS, est d’utiliser différents outils classiques d’analyse (assemblage, mapping, visualisation) et de vous permettre une première prise en main concrète de ce type de données. Nous ne balayerons pas tous les types de données ou toutes les analyses. De manière générale, ne vous contentez pas d’appliquer les instructions données dans le TP, essayez d’explorer les différentes options de chaque outil, en utilisant son aide (généralement présente en bas de page Galaxy de l’outil).

Connectez vous sur l’instance Galaxy Migale. Les login et mot de passe sont les mêmes que ceux de votre poste fixe (vous pouvez également utiliser votre compte Migale). Les données du TP sont accessibles dans le menu Shared Data/Data Libraries/NGS Formation : Un répertoire Reads pour les lectures brutes, un répertoire Refs pour le génome de référence. Nous vous conseillons de créer un historique pour chaque partie du TP afin de séparer vos fichiers de sortie. Dans la mesure du possible, renommez les fichiers de sortie pour leur donner des noms explicites.

2 Logiciels

Nous utiliserons plusieurs programmes pour évaluer la qualité du séquençage, écourter les lectures, assembler, évaluer les assemblages, aligner et visualiser l’alignement. Voici la liste des programmes, accompagnés d’une brève description de leur fonction :

• FASTQC : Outil graphique produisant un rapport synthétique sur la qualité de données FASTQ.
• Sickle : Outil de filtrage et nettoyage des bases mal séquencées, repérées selon des critères de qualité, dans les lectures.
• SPADES : assembleur de données de séquence haut débit.
• QUAST : Calcule un ensemble de mesures sur l’assemblage.
• Bowtie2 : Outil d’alignement des lectures sur un génome de référence.
• IGV : Visualisateur de génome spécialisés dans la représentation des données de mapping.

3 Données, Qualité et Nettoyage

Dans cet TP, nous assemblerons en contigs des séquences d’une souche d’Escherichia coli obtenues à partir d’un séquenceur de paillasse de type MiSeq d’Illumina. Ces données sont mises à disposition par le fabricant. Tous les renseignements sur ce jeu de données sont disponibles sur cette page : https://emea.illumina.com/informatics/sequencing-data-analysis/data-examples.html

Afin d’accélérer les calculs, nous ne travaillerons que sur 10% des données, ce qui permet déjà d’avoir des résultats corrects dans un temps de calcul convenable pour le TP. Votre travail consiste dans un premier temps à effectuer un contrôle qualité sur le jeu de données fourni. Dans un second temps il vous faudra nettoyer les séquences en vous référant aux différentes étapes de nettoyage listées ci-dessous.

3.1 Nettoyage et contrôle qualité

Votre travail consiste dans un premier temps à effectuer un contrôle qualité sur le jeu de données fournis à l’aide du logiciel fastqc [1].

3.1.1 FastQC

• Créez un nouvel historique.
• Importez y les fichier correspondant aux lecture brutes :
  • R1_10.fastq
  • R2_10.fastq
• Lancez le contrôle qualité sur chacun de ces fichiers avec l’outil FastQC présent dans le dossier FASTQ Quality control
  • On peut sélectionner plusieurs dataset avec le bouton “Multiple datasets”.
• Pour la visualisation des résultats, ouvrez la page générée (WebPage). Analysez le rapport produit par FastQC. Le séquençage vous paraît-il de bonne qualité ?
  • Avec certains navigateur, il peut être nécessaire de télécharger la WebPage pour la visualiser correctement.

3.1.2 MultiQC

L’outil multiqc [2] permet d’agréger les résultats de plusieurs outils d’analyses bioinformatiques (dont FastQC) en un seul rapport.

• Regrouper les deux rapports de contrôle qualité FASTQC au format RawData avec l’outil multiQC présent dans le dossier FASTQ Quality control.
  • Si vous avez une erreur, avez vous bien sélectionné le bon outils “Which tool was used generate logs?” ?
• Visualisez (après téléchargement si besoin) et analysez le rapport produit par multiQC. Quelles variations observez-vous entre les deux jeux de données ?

3.1.3 Nettoyage

Le trimming, ou nettoyage des lectures de mauvaise qualité est une étape indispensable avant toute analyse. Cela permet de diminuer le nombre de lectures à traiter et de ne conserver que les parties des lectures de bonne qualité pour l’analyse. Il existe de nombreuses façon de filtrer ces lectures, nous allons utiliser une méthode de trimming “adaptative”, qui analyse les lectures individuellement (ou par paires). Trimmer les lectures brutes en gardant les lectures de taille (après trimming) supérieure à 50 bases, ne comprenant pas de N et dont la qualité moyenne (après trimming) et supérieure à 20, à l’aide de l’outil sickle. Attention, sickle peut être utilisé en mode “single” ou “paired”. Nous tenons à nettoyer les lectures en mode pairé. Sickle doit garder les paires de bonne qualité et mettre dans un fichier à part les lectures de bonnes qualité dont la paire a été filtrée.

• Lancez l’outil Sickle du dossier FASTA/FASTQ en mode pairé avec les paramètres appropriés sur le couple de fichiers bruts de lectures :
  • Paired-end (two separate input file)
    • forward : R1_10.fastq
    • reverse : R2_10.fastq
  • Quality threshold = 20
  • Length threshold = 50
  • Truncate sequences with Ns at first N position = Yes
• Si les fichiers d’entrée ne vous sont pas proposés dans le formulaire c’est que les fichiers fastq ne sont pas automatiquement reconnus comme étant au format fastq-sanger. Vous pouvez, soit utiliser Fastq Groomer pour les convertir, soit si vous êtes sur du format de votre fichier (ce qui est le cas ici), en changer directement le type dans les propriétés de chaque fichier (icône de crayon).
• Renommez les fichiers de sortie :
  • R1Trimmed.fastq
  • R2Trimmed.fastq
  • singletons.fastq
• Relancez FastQC sur les fichiers trimmés
• Comparez les résultats entre fichier trimmés et bruts et observez l’effet du nettoyage.

Questions

• Quelle est le nombre de lectures avant et après nettoyage ?

  Avant : \(250000*2 = 500000\)
  Après : \((236771*2) + 12180 = 485722\)

• Quelle est la longueur moyenne des lectures avant et après nettoyage ?

  Avant : \(151\)
  Après : R1: \(143\), R2: \(134\)

• Sachant que le génome à assembler fait environ 4.6 Mb, déduisez du nombre de bases¹, la profondeur approximative de séquençage avant/après trimming ?

  Avant : \(\frac{250000*2*151}{4600000}=16.41\mathsf{X}\)
  Après : \(\frac{236771*143 + 236771*134}{4600000}=14.26\mathsf{X}\)

4 Assemblage de Novo

4.1 Assemblage avec SPADES

Il existe de nombreux outils dédiés à l’assemblage de novo. SPADES [3] est actuellement l’un des outils les plus utilisé et est considéré à l’état de l’art pour l’assemblage de génomes bactériens. Il est possible de l’utiliser en ligne de commande ou bien avec Galaxy. Spades comprends une phase de correction des reads suivi d’une phase d’assemblage. Il est possible de l’utiliser avec des reads court (comme lors de ce TP) ou en assembleur hybride avec des reads longs.

• Explorer les paramètres et lancer l’assemblage de Spades sur nos données trimmées pairées
• Appliquer les paramètres suivant :
  • Set coverage cutoff option = Auto
  • Select k-mer detection option = Auto
  • Library type : Paired-end (separate input files)
    • R1Trimmed.fastq
    • R2Trimmed.fastq

4.2 Calcul de métriques

Il y a fondamentalement 2 types de métriques calculables sur un assemblage : celles basées sur un génome de référence: couverture, erreurs (ou différences) d’assemblage, … et celles sans génome de référence: taille de l’assemblage, nombre de contigs, N50, …

Rappelons que le N50 correspond à la plus petite taille de contig/scaffold telle que 50% de l’assemblage est contenu dans des contigs/scaffolds plus grands que cette taille.

Le NG50 est la plus petite taille de contig/scaffold telle que la somme des tailles des contigs/scaffolds plus grands que cette taille dépasse 50% de la taille du génome.

Les métriques avec génome de référence sont essentiellement utiles pour l’évaluation des assembleurs. Il est souvent impossible de les calculer précisément quand la référence n’est pas connue. Pour calculer des métriques et générer un rapport complet (incluant les erreurs ou différence d’assemblage si on lui fournit une référence), nous utilisons l’outil quast [4], en prenant comme référence la version Sanger de la souche que nous avons séquencée.

• Importez dans votre historique le génome de référence : ref.fna
• Lancez le calcul de métriques sur l’assemblage (fichier contigs.fa) avec l’outil Quast présent dans la section Assembly en donnant comme référence ref.fna:
  • Contigs/scaffolds file : SPADES ... contigs (fasta)
  • Type of assembly : Genome
  • Use a reference genome? : Yes
    • ref.fna
    • Type of organism : prokaryotes

Questions

• Combien de contigs comprend l’assemblage ?

  223 contigs
  121 contigs (>= 500 bp)

• Quel est le N50 de l’assemblage que vous avez produit ?

  87 059 nct

• Quelle est la taille du plus grand contig ?

  221 645 nct

• Quelle est la taille totale de l’assemblage ?

  4 559 360 nct (>= 500 bp) , pour 4,6 Mb attendus.

• Quel est le NG50 de l’assemblage ?

  73 169 nct

• Quel est le nombre de erreurs/différences par rapport à la référence ?

  80 mismatches, 10 indels et 0 misassemblies.

• Quel fraction du génome de référence est réassemblée ?

  98.185%

• Comparer le %GC de notre assemblage à la référence.

• Ouvrir icarus.html > Contig alignment viewer pour visualiser l’alignement des contigs avec la référence.

5 Mapping contre une référence

5.1 Bowtie2

bowtie2 [5] est un logiciel de mapping de lectures très fréquemment utilisé. Il peut par exemple servir à aligner des reads avec un génome de référence, ou dans notre cas à aligner des reads avec des contigs assemblés. L’alignement se fait en deux étapes : la création d’un index puis l’alignement proprement dit. Là encore, Galaxy permet de lancer facilement ces deux étapes. Pour faciliter les post-traitements, nous demanderons à Bowtie de sortir le résultat d’alignement des lectures pairées trimmées au format BAM. La taille (annoncée) d’insert étant de 300, nous demanderons à bowtie de valider les lectures comme étant pairées si elles ont une taille d’insert comprise entre 200 et 400. Récupérer également les statistiques de mapping.

• Créer un nouvel historique.
• Y copier les fichiers des reads trimmées et la référence.
• Lancer l’alignement avec l’outil bowties2 de la section Mapping pour des fichiers pairés R1Trimmed.fastq et R2Trimmed.fastq sur la référence contigs.fa.
  • Is this single or paired library : paired-end
    • File #1 : R1Trimmed.fastq
    • File #2 : R2Trimmed.fastq
  • Do you want to set paired-end options? : Yes
    • Set the minimum fragment length for valid paired-end alignments : 200
    • Set the maximum fragment length for valid paired-end alignments : 400
  • Will you select a reference genome from your history or use a built-in index? Use a genome from the history and build index : SPADES ... contigs (fasta)
  • Do you want to use presets? No, just use defaults
  • Save the bowtie2 mapping statistics to the history : Yes

• Si vous visualisiez le fichier SAM, pour traduire les drapeaux de SAM (2e colonne), vous pouvez utiliser le site https://broadinstitute.github.io/picard/explain-flags.html.

5.2 Samtools

samtools [6] est un ensemble d’outils permettant de manipuler les fichiers issus d’un alignement. samtools travaille sur des fichiers de type SAM/BAM. Le format BAM est une version binaire du format SAM, qui lui est un format textuel (lisible par les humains). samtools propose plusieurs outils en ligne de commande pour transformer les fichiers SAM en BAM, faire des statistiques, visualiser un alignement ou encore générer la distance entre la lecture et la séquence de référence sur lequel il est aligné. Beaucoup d’outils prennent du BAM en entrée. On peut passer de SAM à BAM (et inversement) sans difficultés.

Un certain nombre d’outils de samtools sont interfacés dans Galaxy. Nous allons en utiliser quelques uns :

5.2.1 Création d’un fichier BAM

Lorsque l’on utilise les outils d’alignement en ligne de commande, il est important de convertir le SAM en BAM puis de le trier et de l’indexer. Ici, ces étapes sont gérées par Galaxy et transparentes pour nous.

5.2.2 Statistiques sur les alignements

• Utiliser l’outil Samtools flagstats qui prend en entrée le fichier BAM trié pour calculer des statistiques sur l’alignement.

Questions

• Quel est le nombre de lectures alignées pairées (properly paired) ?

  466386 properly paired (98.49%)

• Utiliser l’outil Samtools idxstats qui prend lui aussi en entrée le fichier BAM trié pour calculer des statistiques de profondeur?

Questions

• Combien de reads sont alignés sur le premier contigs ?

  22520 reads

5.3 Visualisation

Il existe des outils de visualisation directement hébergé par Galaxy, cependant ils ne sont pas très stables. Nous allons donc privilégier IGV (Integrative Genomics Viewer).

• Télécharger les fichiers
  • le fichier de reads alignés et triés BAM
  • l’index BAI
  • la séquence de référence (c’est à dire des contigs) contigs.fa
• Lancer IGV sur votre poste personnel
• Charger la référence : Genomes > Load Genome from File > contigs.fa
• Charger les reads : File > Load from File > *.bam (pour rappel, il faut que l’index soit disponible dans le même dossier avec le même nom et l’extension .bai)

6 Et sur ses données, … comment fait on ?

Il faut bien suivre les étapes montrées dans ce TP, notamment le contrôle qualité pour s’assurer qu’il n’y a pas de biais inattendu dans ses données. Concernant l’assemblage, nous vous conseillons d’utiliser un assembleur à l’état de l’art , comme SPADES. Vous pouvez utiliser des “optimiseurs” d’assemblage comme Unicycler ou Shovill qui cherchent à optimiser les paramètres d’assemblage, circulariser vos génomes et améliorer ainsi la qualité globale de vos assemblages.

7 Outils, manuels et sites d’interêt

7.1 Outils

• FastQC : https://www.bioinformatics.babraham.ac.uk/projects/fastqc/
• Sickle : https://github.com/najoshi/sickle
• prinseq : http://prinseq.sourceforge.net/
• SPAdes : https://cab.spbu.ru/software/spades/
• QUAST : https://quast.bioinf.spbau.ru
• Bowtie2 : https://bowtie-bio.sourceforge.net/bowtie2/index.shtml
• SAMtools : http://www.htslib.org/
• IGV : https://software.broadinstitute.org/software/igv/

7.2 Sites

• https://biostars.org : Forum de questions/réponses généralistes bioinformatique. Très actif côté analyse des données NGS.
• https://seqanswers.com : Forum plus porté sur la partie biologique des NGS

Les références

1. Andrews S. FastQC A Quality Control tool for High Throughput Sequence Data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/. 2010. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.

2. Ewels P, Magnusson M, Lundin S, Käller M. MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics. 2016;32:3047‑8.

3. Bankevich A, Nurk S, Antipov D, Gurevich AA, Dvorkin M, Kulikov AS, et al. SPAdes: A New Genome Assembly Algorithm and Its Applications to Single-Cell Sequencing. Journal of Computational Biology. 2012;19:455‑77. doi:10.1089/cmb.2012.0021.

4. Gurevich A, Saveliev V, Vyahhi N, Tesler G. QUAST: quality assessment tool for genome assemblies. Bioinformatics. 2013;29:1072‑5. doi:10.1093/bioinformatics/btt086.

5. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nature methods. 2012;9:357.

6. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078‑9.

Notes de bas de page

1. Galaxy ne permet pas de calculer de façon simple le nombre de bases dans un fichier fastq. Contentez vous donc d’une approximation en utilisant le nombre de reads et la longueur moyenne. Pour ceux qui voudraient le chiffre exact, vous pouvez le calculer en utilisant les outils FASTQ to Tabular suivie de Cut puis Line/Word/Character count.