Alignement de séquences avec BLAST
Cet article a pour but de vous montrer une application pratique de BLAST, le fameux programme d'alignement de séquences détenant un record de citations, avec certains problèmes qu'on peut rencontrer et ce qu'on peut tirer de son résultat. BLAST a au moins deux usages typiques en génomique :
- Trouver les occurrences similaires à une séquence de bases ou d'acides aminés donnée ("query") dans une séquence de référence, par exemple un génome entier.
- Comparer deux séquences supposées semblables (par exemple des gènes orthologues ou leurs protéines) pour mesurer leur similarité et repérer ce qui est conservé/muté.
Pour les initiés, le premier point est à ne pas confondre avec ce qu'on utilise en séquençage à haut débit pour aligner des millions de petites séquences. BLAST recherche des séquences similaires, et pas l'emplacement de la séquence exacte (à erreurs de séquençage près). Il alignera une ou plusieurs parties de la séquence, éventuellement séparées par des "trous" ("gaps"). Le principal critère pour qu'il reporte un alignement est qu'il arrive à trouver au moins un "mot" de taille 11 (par défaut, on peut le changer) présent exactement dans la référence, puis il étend l'alignement autour de ce mot.
BLAST possède une interface web facile d'utilisation, même s'il faudra comprendre l'algorithme de Smith-Waterman pour être en mesure d'ajuster ses paramètres avancés. Mais nous allons l'utiliser en ligne de commande pour plus d'efficacité et de stabilité, et voir comment retirer du résultat les infos qui nous intéressent.
Un exemple pratique
Je vous propose l'application suivante, tirée d'un cas réel récent :
- On recherche près des télomères (à moins de 500Kb) des domaines appelés CRISPR, qui serviront de cible pour insérer artificiellement des fragments d'ADN dans l'organisme. En voilà un exemple : GCGGGTCCTGTAGTGCCGGG.
- Les régions télomériques de l'assemblage de référence habituel (de NCBI) sont incorrectes, et le "client" possède sa propre séquence de référence, plus précise.
- Le client veut trouver, parmi une douzaine de CRISPR, lequel à la fois se trouve sur un maximum de régions télomériques, et est absent du reste du génome. C'est-à-dire, lequel est le plus spécifique aux régions d'intérêt (pour éviter d'insérer de l'ADN codant n'importe où ailleurs).
Partons du principe que vous avez installé BLAST. Commençons par préparer la "query" : les séquences CRISPR à rechercher. Pour cela on crée un fichier "crispr.fasta", de ce type :
|
1 2 3 4 5 6 |
>seq1 GTGCCCCGGCGCCACGANGG >seq2 GTGCCCGGCTGCAAGCANGG >seq3 CGCCCCCTCGTGGCGCCNGG |
Notez le nucléotide "N" (="n'importe quelle base"). En principe, BLAST devrait gérer les consensus IUPAC, mais en fait ça peut être très limitant en pratique (voir ici, 7.2).
Ensuite, comme nous fournissons nous-mêmes la référence "reference.fasta" (pour l'exemple, on pourrait prendre le chromosome 20 de l'humain : chr20.fa.gz, à décompresser et renommer), il faut créer l'index correspondant avec la commande "formatdb" en ligne de commande :
|
1 |
formatdb -p F -i reference.fasta -o T |
ce qui créera plein de petits fichiers mystérieux, que vous préférerez peut-être grouper dans un dossier propre. Le "-p F" indique qu'on aligne des nucléotides ("-p T" pour les protéines). Le "-o T" est moins clair, lié au format des headers du fasta, mais conseillé tant que ça marche.
Enfin on peut lancer BLAST en lui donnant la query et la référence, comme ça :
|
1 |
blastall -p blastn -d reference.fasta -i crispr.fasta -o alignments.txt |
"blastall" est le programme générique, puis on peut choisir le type d'alignement : "blastn"/"nc_megablast"/"megablast" (nucléotides), "blastp" (acides aminés), etc. Nous utilisons le standard "blastn", très permissif, mais si on cherchait des séquences très proches de la nôtre, on préférerait "megablast" ou "nc_megablast" (nc_ pour non-continuous) - voir ici comment choisir le bon.
Le résultat, "alignments.txt", peut être assez lourd. Pour chaque séquence (par ex. chromosome) de la référence, il donne tous les alignements trouvés de chaque séquence de crispr.fasta. Un extrait correspondant à la query "seq3" (en utilisant le chromosome 20 de l'humain - hg19 - comme référence) :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
[...] >chr20 Length = 63025520 Score = 30.2 bits (15), Expect = 7.4 Identities = 15/15 (100%) Strand = Plus / Minus Query: 1 cgccccctcgtggcg 15 ||||||||||||||| Sbjct: 1946483 cgccccctcgtggcg 1946469 Score = 30.2 bits (15), Expect = 7.4 Identities = 15/15 (100%) Strand = Plus / Plus Query: 3 ccccctcgtggcgcc 17 ||||||||||||||| Sbjct: 34330497 ccccctcgtggcgcc 34330511 [...] |
On voit que "seq3" a trouvé deux homologues sur le chromosome 20 : les bases 1-15 en position 1946483-1946469, et les bases 3-17 en 34330497-34330511. Un score d'enrichissement est donné, ainsi que le nombre de matches sans erreur ("Identities"). Tout en haut est indiquée la taille du chromosome ; on voit ainsi que ces deux alignements particuliers sont loin des télomères, donc pas du tout ce qui nous intéresse.
Mais quand il y a des milliers de lignes comme celles-ci, on se dit qu'il faut un moyen plus systématique d'examiner les résultats. De toute évidence, le parser soi-même est un cauchemar ; heureusement, des outils sont déjà disponibles. J'ai choisi Biopython, la librairie bioinfo du langage Python. Voilà comment examiner le contenu du fichier :
|
1 2 3 4 5 6 |
from Bio import SearchIO records = SearchIO.parse("alignments.txt", "blast-text") # iterateur for record in records: # record est un objet QueryResult for hit in record: # hit est un objet Hit for hsp in hit: # hsp est un objet HSP print hsp |
Les objets QueryResult, Hit et HSP contiennent et rendent accessible toute l'info du fichier à travers leurs attributs. Un QueryResult concerne une query, en particulier :
|
1 2 3 4 5 6 |
record.hit_keys # les references touchees par cette query record.hsps # raccourci pour iterer sur les HSPs record.id # nom de la query (comme "seq1") record.program # programme utilise (comme "blastn") record.seq_len # longueur de la query (20) record.target # nom de la sequence de reference |
Un Hit est l'ensemble des alignements de la query sur une certaine référence, en particulier :
|
1 2 3 4 |
hit.hsps # liste des HSPs hit.id # nom de la reference touchee (comme "chr1") hit.query_id # nom de la query hit.seq_len # longueur de la reference touchee |
Un HSP est un alignement particulier. Il a de nombreux attributs qui parlent d'eux-mêmes pour la plupart :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
>>> dir(hsp) [..., 'aln', 'aln_all', 'aln_annotation', 'aln_annotation_all', 'aln_span', 'alphabet', 'bitscore', 'bitscore_raw', 'evalue', 'fragment', 'fragments', 'gap_num', 'hit', 'hit_all', 'hit_description', 'hit_end', 'hit_end_all', 'hit_features', 'hit_features_all', 'hit_frame', 'hit_frame_all', 'hit_id', 'hit_inter_ranges', 'hit_inter_spans', 'hit_range', 'hit_range_all', 'hit_span', 'hit_span_all', 'hit_start', 'hit_start_all', 'hit_strand', 'hit_strand_all', 'ident_num', 'is_fragmented', 'pos_num', 'query', 'query_all', 'query_description', 'query_end', 'query_end_all', 'query_features', 'query_features_all', 'query_frame', 'query_frame_all', 'query_id', 'query_inter_ranges', 'query_inter_spans', 'query_range', 'query_range_all', 'query_span', 'query_span_all', 'query_start', 'query_start_all', 'query_strand', 'query_strand_all'] |
Ce qui nous permet de compter par exemple combien de fois "seq1" s'aligne à plus de 500Kb des télomères :
|
1 2 3 4 5 6 7 8 9 10 11 |
from Bio import SearchIO records = SearchIO.parse("alignments.txt", "blast-text") record = records.next() # la premiere query intergenomic = 0 telomeric = 0 for hit in record: for hsp in hit: if hsp.hit_range[1] > 500000 and hsp.hit_range[0] < hit.seq_len-500000: intergenomic += 1 else: telomeric += 1 |
Résultat :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
>>> print "Nombre de hits intergenomiques:", intergenomic 0 >>> print "Nombre de hits telomeriques", telomeric 1 >>> print hsp # le seul alignement, encore dans cette variable Query: seq1 Hit: chr20 Query range: [1:14] (1) Hit range: [62917482:62917495] (1) Quick stats: evalue 4.7; bitscore 26.30 Fragments: 1 (13 columns) Query - tgccccggcgcca ||||||||||||| Hit - tgccccggcgcca |
Sur ce chromosome au moins, tout va bien ! "seq1" n'est présente qu'une fois près du télomère, et nulle part ailleurs, même si seulement 13 nucléotides s'alignent.
A vous ensuite d'adapter le script à vos besoins particuliers.
Remarques :
- On peut aussi demander à BLAST une sortie en format XML.
- On peut aussi lancer BLAST directement depuis Biopython !
Merci à Sylvain P., Clem_, ook4mi et Yoann M. pour leur relecture.
Jacques D
mai 9, 2014 à 2:05
J’aime beaucoup cet article pour son esprit pratique. Le Blast et le “parsage” de ses resultats sont des elements recurrent en bioinfo. Il est donc bon de donner des fiches pratiques pour tout nouveau utilisateur (bioinformaticien et autres) ou rafraichir nos memoires trop souvent faillibles.
Je souhaiterai contribuer au sujet:
- D’abord une information utile pour les cas des courtes sequences comme celle de 20 nucleotides du cas pratique decrit ici. Le “Word size” par defaut, par exemple sur le site du ncbi, est de 28. Ce qui signifie qu’il sera impossible de trouver des sequences qui s’aligne sur moins de 28 nucleotides… Mais ils ont pense a tout: “Instead, the nucleotide and protein blast programs automatically check for short queries and adjust the search parameters accordingly. This adjustment occurs when the query, either nucleotide or amino acid, is of length 30 or less.” Donc dans ce cas on aura bien un resultat pour cette sequence. Apparement la version de blast que tu utilise ici en local arbore le meme comportement.
La version “Nucleotide-Nucleotide BLAST 2.2.28+” que j’utilise sur mon Mac semble avoir un “Word size” par defaut a 56 et ne s’auto-ajuste pas. Du coup, dans un cas comme le mien, pour faire ton etude il faut penser a reduire le parametre (max 15 … min 4) pour obtenir un resultat.
- Ensuite j’aime beaucoup utiliser le bash pour se genre d’etude. Je vous presente une maniere de proceder pour avoir les sequences pres des telomeres qui sont absentes des regions intergenomiques:
1) blastn -db reference.fa -query crisp.fa -outfmt '6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send slen' -word_size 4 > resu.txt
Cette commande ecrit les resultats du blast au format tabule dans le fichier resu.txt. (Le format tabule se choisit avec l’option 6. Vous pouvez vous-meme choisir d’autres info ou changer leurs ordres. Voir: blastn -help”)
2) seqs=$(awk '{if (($9 $11-500000) && ($10 $11-500000)) {telomeric[$1]="x"} else {intergenomic[$1]="x"} } END { for (i in telomeric){ if (intergenomic [i] == "") print i}}' resu.txt);
for i in $seqs; do
awk -v test=$i '{if ($1 == test ) print $0}' resu.txt
done
La commande awk peut sembler complexe (pour rappel, $0=toute la ligne courante, $1=colonne 1 de la ligne courante, $2= colonne 2 de la ligne courante, etc. Et bien sur il parcourt les lignes une par une de haut en bas).
Pour filtrer les resultats de blasts awk est tres utile. Par exemple, dans mon resultat mon pourcentage d’identite (pident) se retrouve dans la colonne 3. Pour afficher seulement les resultats de blast qui ont un pident > 95 il suffit de faire :
awk '{if($3>95) print $0}' resu.txt
on peut ainsi rajouter autant de test que necessaire …
P.S : Je sais pas si ca existe deja sur le site, mais il serait peut-etre bon d’ajouter aussi l’approche en BioPerl…
P.S2 : Desole j’ai pas d’accent sur mon clavier …
Julien Delafontaine
mai 9, 2014 à 2:42
Merci pour ce commentaire !
D'apres ce que j'ai lu, c'est Megablast qui a un word size de 28, et blastn est a 11 :
"The most important reason that blastn is more sensitive than MEGABLAST is that it uses a shorter default word size (11)". (http://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=ProgSelectionGuide#tab31, 4.2).
Par contre c'est vrai qu'il faut bien penser a l'ajuster car c'est le parametre le plus important, par exemple un blastn avec des sequences de longueur >100.
Merci aussi pour la version bash. Je ne connais pas assez Perl mais je pourrais ecrire une version R si quelqu'un s'y interesse.
Jacques D
mai 9, 2014 à 3:42
C'est une bonne idée une version R !
Autant pour moi, pour le "Word Size" tu as raison. Par défaut il te mette sur megablast quand tu choisis blastn sur NCBI ... d'ou ma confusion.
Par contre, il est intéressant de noter que megablast a ce comportement adaptatif en fonction de la longueur de la séquence, dont je parlais précédemment (http://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=FAQ).
Par contre je comprends toujours pas pourquoi ma version (telecharge sur le NCBI) a un Word Size de 56 ... le côte obscur de la bioinfo surement !
P.S: les accents ont été ajoute par auto-completion ... pas terrible.
Mousset Carla
mars 19, 2016 à 10:49
Je suis étudiante en Master 2 en Ecologie et j'aimerai s'il vous plait savoir comment faire pour choisir la bonne séquence similaire résultat après avoir lancé une requête sur BLAST en ligne.
Pour être plus précise, j'ai des résultats de séquençage pour lesquels j'ai des séquences de nucléotide identifiées à l'ordre par exemples, dont plusieurs propositions d'espèces ont été faites comme résultat de l'identification. Pour ces cas lorsque je blast la séquence j'ai effectivement toutes ces propositions d'espèces toutefois, je ne sais pas comment choisir la bonne séquence.
Merci d'avance pour votre réponse
Julien Delafontaine
avril 28, 2016 à 5:48
Bonjour, malheureusement je ne me souviens plus du tout de tout ça, mais sur la version en ligne, je me souviens que les résultats apparaissent triés de haut en bas du meilleur au moins bon, et ceci est possible parce que chaque alignement a un score, comme dans mon exemple où on peut lire " Score = 30.2 bits (15), Expect = 7.4".
coraline
mai 13, 2016 à 9:30
Bonjour,
Je suis étudiante en l2 génie biologie option bio-informatique, je me demandais s'il était possible de lancer en local une commande permettant d'obtenir la taille d'un fragment entre deux amorces.
Je m'explique, je connais mes deux amorces, reward et reverse et j'ai téléchargé la séquence de la bactérie qui m'interesse, je voudrais donc pouvoir faire matcher mes deux amorces avec cette séquence_ref et obtenir la taille de la séquence entre les deux amorces (ici ça me permettrai d'otenir la taille de l'ARNr S16) dans un fichier result.
Est-ce possible et si oui pourriez vous m'aider pour la commande ?
Merci d'avance pour votre réponse,
Cordialement,
Coraline.
rk
avril 2, 2020 à 7:43
Bonjour ya t'il un moyen de bous contacter ? C'est au sujet de la génétique..
Yoann M.
avril 3, 2020 à 9:23
Bonjour rk,
Pour nous contacter nous vous invitons à nous envoyer un mail à l'adresse qui figure en bas de page d'un article : "Pour tout renseignement supplémentaire : admin AT bioinfo-fr.net"
Ismael
août 2, 2020 à 1:12
Bonjours
Je suis étudiant en master biochimie, je voudrai savoir comment les résultats de blast sont-ils organisés?