Sciences du vivant

Mikaël Salsonmikael.salson@univ-lille.fr

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.

Thème 1 – algorithmique pour l'analyse de données de séquençage à haut débit

Le séquençage à haut débit de l'ADN, comment ça fonctionne ?

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

Ressources vidéos

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).

Pourquoi vouloir retrouver les reads dans un génome ?

Why do We Map Reads?

Un annuaire pour génome !

Cours Algorithms for DNA sequencing sur Coursera : Indexing and the k-mer index

Introduction au principe de l'alignement

The Alignment Game and the Longest Common Subsequence Problem

Programmation dynamique et calcul d'alignement

Dynamic Programming and Backtracking Pointers

Alignement de séquences global et local

From Global to Local Alignment

Des diapos de cours

Sur l'alignement de séquences et la programmation dynamique : 1, 2, 3

Un article complémentaire

Les structures d'indexation « classiques » et les vecteurs de bits et wavelet tree

Stockage efficace et « simple » de données redondantes

Un article sur un index utilisant les block trees.

Thème 2 – Rechercher dans des milliers de jeux de séquençage

Population BWT
Quantité de données stockée dans SRA
Using reference-free compressed data structures to analyse sequencing reads from thousands of human genomes
Filtres de Bloom
Explication et exemples interactifs des filtres de Bloom
Vidéo sur une analyse des filtres de Bloom
Page Wikipedia assez complète sur les filtres de Bloom
Indexer avec des filtres de Bloom
Sequence Bloom Trees
Split Sequence Bloom Trees
Le SSBT expliqué sur le blog de Carl Kingsford (un des auteurs)
Le AllSome Sequence Bloom Trees

Thème 3 – Recherche dans des données bruitées

Activité 1 – Lecture d'articles

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…).

Comment indexer et interroger des centaines de milliers de données de séquençage ?

Sub-linear Sequence Search via a Repeated And Merged Bloom Filter (RAMBO)

Comment compter les occurrences d'une séquence dans des centaines de milliers de données de séquençage ?

REINDEER: efficient indexing of k-mer presence and abundance in sequencing datasets

Comment classer de longues séquences en fonction de leur espèce d'appartenance ?

Strain-level metagenomic assignment and compositional estimation for long reads with MetaMaps

Comment identifier des variants génétiques dans des milliers de génomes

Comment indexer des milliers de génomes ?

Haplotype-aware graph indexes

Comment stocker un ensemble de manière approximative ?

Xor Filters: Faster and Smaller Than Bloom and Cuckoo Filters

Comment stocker efficacement des vecteurs de bit compressés ?

Tree-Encoded Bitmaps"

Activité 2 – mise en pratique

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.

Se familiariser avec les données

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.

Spécification

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.

Indications pour l'indexation

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.

  1. Solution naïve : parcours de L
  2. Précalcul du rank pour des blocs de taille fixe : cela évite de parcourir totalement L, il suffit de connaître le rank au début du bloc (précalculé) et ensuite de parcourir un bloc jusqu'à la position d'intérêt. C'est assez simple à implémenter et plus efficace que la solution naïve
  3. Utilisation d'un wavelet tree, soit en implémentant soi-même les vecteurs de bits soit en utilisant une bibliothèque externe pour cela

Rendu

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.

Jeux de données

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.

Évaluations des années précédentes