![alignement_blast](https://i0.wp.com/bioinfo-fr.net/wp-content/uploads/2014/04/Screen-Shot-2014-04-15-at-10.51.35-AM.png?resize=433%2C338&ssl=1)
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 7 |
> ;seq1<br> GTGCCCCGGCGCCACGANGG<br> > ;seq2<br> GTGCCCGGCTGCAAGCANGG<br> > ;seq3<br> CGCCCCCTCGTGGCGCCNGG<br> Notez le nucléotide "N" (="n'importe quelle base"). En principe, BLAST devrait gérer les  ; consensus <a href="http://www.bioinformatics.org/sms/iupac.html" target="_blank" rel="noopener">IUPAC</a>, mais en fait ça peut être très limitant en pratique (voir <a href="http://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=ProgSelectionGuide" target="_blank" rel="noopener">ici</a>, 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 2 |
formatdb -p F -i reference.fasta -o T<br> 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 <em>headers</em> 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 2 |
blastall -p blastn -d reference.fasta -i crispr.fasta -o alignments.txt<br> "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 <em>non-continuous</em>) - voir <a href="http://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=ProgSelectionGuide" target="_blank" rel="noopener">ici</a> 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 |
[...] |
>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 7 |
from Bio import SearchIO<br> records = SearchIO.parse("alignments.txt", "blast-text") # iterateur<br> for record in records : # record est un objet QueryResult<br> for hit in record : # hit est un objet Hit<br> for hsp in hit : # hsp est un objet HSP<br> print hsp<br> Les objets QueryResult, Hit et HSP contiennent et rendent accessible toute l'info du fichier à travers leurs attributs. Un QueryResult concerne une <em>query</em>, en particulier : |
1 2 3 4 5 6 7 |
record.hit_keys # les references touchees par cette query<br> record.hsps # raccourci pour iterer sur les HSPs<br> record.id # nom de la query (comme "seq1")<br> record.program # programme utilise (comme "blastn")<br> record.seq_len # longueur de la query (20)<br> record.target # nom de la sequence de reference<br> Un Hit est l'ensemble des alignements de la <em>query</em> sur une certaine référence, en particulier : |
1 2 3 4 5 |
hit.hsps # liste des HSPs<br> hit.id # nom de la reference touchee (comme "chr1")<br> hit.query_id # nom de la query<br> hit.seq_len # longueur de la reference touchee<br> 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 16 |
> ;> ;> ; dir(hsp)<br> [..., 'aln', 'aln_all', 'aln_annotation', 'aln_annotation_all',<br> 'aln_span', 'alphabet', 'bitscore', 'bitscore_raw', 'evalue',<br> 'fragment', 'fragments', 'gap_num', 'hit', 'hit_all',<br> 'hit_description', 'hit_end', 'hit_end_all', 'hit_features',<br> 'hit_features_all', 'hit_frame', 'hit_frame_all', 'hit_id',<br> 'hit_inter_ranges', 'hit_inter_spans', 'hit_range', 'hit_range_all',<br> 'hit_span', 'hit_span_all', 'hit_start', 'hit_start_all',<br> 'hit_strand', 'hit_strand_all', 'ident_num', 'is_fragmented',<br> 'pos_num', 'query', 'query_all', 'query_description', 'query_end',<br> 'query_end_all', 'query_features', 'query_features_all',<br> 'query_frame', 'query_frame_all', 'query_id', 'query_inter_ranges',<br> 'query_inter_spans', 'query_range', 'query_range_all',<br> 'query_span', 'query_span_all', 'query_start', 'query_start_all',<br> 'query_strand', 'query_strand_all']<br> 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 12 |
from Bio import SearchIO<br> records = SearchIO.parse("alignments.txt", "blast-text")<br> record = records.next() # la premiere query<br> intergenomic = 0<br> telomeric = 0<br> for hit in record :<br> for hsp in hit :<br> if hsp.hit_range[1] > ; 500000 and hsp.hit_range[0] < ; hit.seq_len-500000 :<br> intergenomic += 1<br> else :<br> telomeric += 1<br> Résultat : |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
> ;> ;> ; print "Nombre de hits intergenomiques:", intergenomic<br> 0<br> > ;> ;> ; print "Nombre de hits telomeriques", telomeric<br> 1<br> > ;> ;> ; print hsp # le seul alignement, encore dans cette variable<br> Query : seq1<br> Hit : chr20<br> Query range : [1 :14] (1)<br> Hit range : [62917482 :62917495] (1)<br> Quick stats : evalue 4.7 ; bitscore 26.30<br> Fragments : 1 (13 columns)<br> Query - tgccccggcgcca<br> |||||||||||||<br> Hit - tgccccggcgcca<br> 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.
Laisser un commentaire