L'objectif de ces six séances est de découvrir, comprendre et appréhender les structures de données et méthodes algorithmiques couramment utilisées en bioinformatique des séquences, en particulier pour l'analyse de millions (ou milliards, ou …) de séquences génomiques.
La vidéo présente la première et la deuxième générations de séquençage de l'ADN. Une troisième génération est en cours d'émergence.
Next Generation Sequencing (NGS) - An Introduction
On appelle read (ou lecture en bon français) les séquences produites par un séquenceur à haut débit. Elles sont typiquement composées de 100 à 1000 nucléotides (voire 10000 pour les plus récentes technologies).
Cours Algorithms for DNA sequencing sur Coursera : Indexing and the k-mer index
The Alignment Game and the Longest Common Subsequence Problem
Dynamic Programming and Backtracking Pointers
From Global to Local Alignment
Sur l'alignement de séquences et la programmation dynamique : 1, 2, 3
Pour cette séance je vous demande de lire un article et d'en faire un résumé oral. Dans le résumé prenez garde à bien présenter les méthodes ou concepts introduits dans l'article mais également le contexte : pourquoi s'intéresse-t-on à cela.
Si des expériences sont faites, pour évaluer les avantages et inconvénients d'un logiciel, identifiez les qualités et les défauts (qui ne sont généralement pas vraiment mis en avant…) du logiciel présenté.
Voici les articles proposés à la lecture. Certains articles ne seront peut-être pas en accès-libre. Dans ce cas-là dirigez-vous vers le site Sci-Hub (à n'utiliser évidemment que si vous avez une copie de l'article dans votre grenier, mais que vous ne l'avez pas avec vous…).
Le projet consiste à indexer tous les génomes du SARS-CoV-2 connus et à analyser des jeux de données de séquençage pour savoir quel génome du virus s'y trouvent.
Il faut donc prévoir une indexation de dizaines de milliers de petits génomes. Le génome du SARS-CoV-2 fait environ 30 000 nucléotides. Les séquences à indexer sont téléchargeables ici au format FASTA.
L'index créé doit permettre de savoir si telle ou telle séquence s'y trouve et de connaître son nombre d'occurrences. On pourra par exemple chercher des MEM (Maximal Exact Matches), c'est-à-dire la séquence la plus longue qu'on retrouve de manière exacte dans l'index. Si on en trouve qu'une seule occurrence on peut considérer qu'on a trouvé un seul génome qui correspond. À l'inverse tant qu'il existe plusieurs occurrences cela signifie que la séquence recherchée n'est pas suffisamment grande pour déterminer de quel souche du virus elle provient.
Les séquences avec lesquelles nous interrogerons l'index seront des reads de séquençage à haut-débit. Ces expériences de séquençage à haut débit peuvent avoir diverses origines : soit des séquençages destinés à identifier la séquence du virus lui-même, soit des séquençages de patients dans lesquels on peut retrouver des traces du virus.
Chaque read est beaucoup plus court que le génome du virus lui-même. On ne s'attend donc qu'à en retrouver une petite partie et tous les reads ne vont pas correspondre à une partie du virus (parfois ce sera peut-être une infime minorité de reads qui correspondront… voire aucun).
Une manière de procéder est de rechercher les MEMs dans les reads et de vérifier s'ils n'apparaissent qu'une fois dans l'index ou ont plusieurs occurrences.
La sortie sera au format TSV et aura la forme suivante :
identifiant1 nombre_d_occurrences identifiant2 nombre_d_occurrences
Pour chaque génome identifié dans les reads en entrée, elle donnera le nombre de fois où il a été vu dans les reads. Dans le cas où il n'a pas été possible d'isoler un MEM unique (donc dans le cas où le MEM aurait plusieurs occurrences), on considérera que chacun des N génomes dans lequel le MEM a été identifié contribue pour 1/N au nombre d'occurrences. On pourra donc avoir des nombres d'occurrences qui ne sont pas entiers.
Votre programme prendra une option qui, si elle est spécifiée, ne considérera que les MEM uniques.
Commencez par analyser les génomes fournis. Il s'agit de tous les génomes récupérés sur le site du NCBI. Or il est tout à fait possible, et même assez probable, que plusieurs génomes soient totalement identiques.
Commencez par identifier les génomes qui sont identiques et à partir de cela créez un nouveau fichier qui ne contienne que des séquences différentes. Lorsque plusieurs séquences sont identiques, attribuez leur un nouvel identifiant et précisez la ou les zones géographiques (comme dans le fichier d'origine). Ainsi une entête d'une séquence après unification pourrait ressembler à :
>UNI00001 | Severe acute respiratory syndrome coronavirus 2 isolate SARS-CoV-2/human | Belgium, Germany
On a donc attribué un identifiant arbitraire à la séquence UNI00001 et on a donné les localisations où ce génome a été observé.
Dans un fichier au format TSV, gardez la trace des unifications que vous avez faites, sous la forme suivante :
UNI00001 MT745000.1 MT745020.1 MT745050.1
Ce fichier servira à indiquer quelles séquences ont été rassemblées sous notre identifiant unique.
Vous produirez un programme index qui prend en argument le fichier FASTA contenant les génomes uniques ainsi que le nom de base (sans extension) de l'index à produire (vous pouvez créer plusieurs fichiers si nécessaire).
Vous produirez un programme query qui prend en paramètre un nom de base d'un index, un fichier FASTQ contenant des reads et un nom de fichier TSV dans lequel seront écrits les résultats (comme indiqué plus haut pour connaître le nombre d'occurrences de chaque génome).
Chacun de ces deux programmes peuvent prendre des options servant à modifier leur comportement. Mais ces options doivent être… optionnelles.
L'indexation devra passer par l'utilisation d'une transformée de Burrows-Wheeler ou d'une table des suffixes. Il faut faire très attention à l'étape de construction afin que celle-ci ne nécessite pas un espace quadratique dans la taille des séquences en entrée.
Vous n'utiliserez pas de bibliothèque externe pour calculer et utiliser la table des suffixes ou la transformée de Burrows-Wheeler.
Dans le cas de l'utilisation d'une transformée de Burrows-Wheeler, différents choix peuvent être faits pour le calcul de la fonction rank sur L.
Votre tavail est à rendre avant le 04 novembre à 20h à : mikael.salson@univ-lille.fr
Vous devez rendre votre code source (évidemment). Un make dans le répertoire de base doit compiler votre programme (s'il y a lieu, sinon il ne fait rien). Un make example doit lancer l'indexation sur un petit exemple puis lancer une requête dessus.
Vous devez également rendre un fichier README (format texte, Markdown ou Org) qui détaille :
Vous ne rendrez surtout pas les jeux de données utilisés.
Votre note dépendra du fonctionnement de votre programme, de son exactitude, de son efficacité, des choix effectués (une méthode plus avancée rapportant plus de points) et de la qualité du rendu.
Les jeux de reads que vous analyserez sont les fichiers au format fastq.gz accessibles ici.
Les noms des fichiers qui finissent par _1 ou _2 correspondent à des séquençages paired-end. Vous pouvez concaténer les deux fichiers.