BLAST en pratique

alignement_blast
Ali­gne­ment de séquences avec BLAST

Cet article a pour but de vous mon­trer une appli­ca­tion pra­tique de BLAST, le fameux pro­gramme d'ali­gne­ment de séquences déte­nant un record de cita­tions, avec cer­tains pro­blèmes qu'on peut ren­con­trer et ce qu'on peut tirer de son résul­tat. BLAST a au moins deux usages typiques en géno­mique :

  1. Trou­ver les occur­rences simi­laires à une séquence de bases ou d'acides ami­nés don­née ("que­ry") dans une séquence de réfé­rence, par exemple un génome entier.
  2. Com­pa­rer deux séquences sup­po­sées sem­blables (par exemple des gènes ortho­logues ou leurs pro­téines) pour mesu­rer leur simi­la­ri­té et repé­rer ce qui est conservé/​muté.

 

Pour les ini­tiés, le pre­mier point est à ne pas confondre avec ce qu'on uti­lise en séquen­çage à haut débit pour ali­gner des mil­lions de petites séquences. BLAST recherche des séquences simi­laires, et pas l'emplacement de la séquence exacte (à erreurs de séquen­çage près). Il ali­gne­ra une ou plu­sieurs par­ties de la séquence, éven­tuel­le­ment sépa­rées par des "trous" ("gaps"). Le prin­ci­pal cri­tère pour qu'il reporte un ali­gne­ment est qu'il arrive à trou­ver au moins un "mot" de taille 11 (par défaut, on peut le chan­ger) pré­sent exac­te­ment dans la réfé­rence, puis il étend l'alignement autour de ce mot.

BLAST pos­sède une inter­face web facile d'utilisation, même s'il fau­dra com­prendre l'algo­rithme de Smith-Water­man pour être en mesure d'ajuster ses para­mètres avan­cés. Mais nous allons l'utiliser en ligne de com­mande pour plus d'efficacité et de sta­bi­li­té, et voir com­ment reti­rer du résul­tat les infos qui nous inté­ressent.

Un exemple pratique

Je vous pro­pose l'application sui­vante, tirée d'un cas réel récent :

  • On recherche près des télo­mères (à moins de 500Kb) des domaines appe­lés CRISPR, qui ser­vi­ront de cible pour insé­rer arti­fi­ciel­le­ment des frag­ments d'ADN dans l'organisme. En voi­là un exemple : GCGGGTCCTGTAGTGCCGGG.
  • Les régions télo­mé­riques de l'assemblage de réfé­rence habi­tuel (de NCBI) sont incor­rectes, et le "client" pos­sède sa propre séquence de réfé­rence, plus pré­cise.
  • Le client veut trou­ver, par­mi une dou­zaine de CRISPR, lequel à la fois se trouve sur un maxi­mum de régions télo­mé­riques, et est absent du reste du génome. C'est-à-dire, lequel est le plus spé­ci­fique aux régions d'intérêt (pour évi­ter d'insérer de l'ADN codant n'importe où ailleurs).

 

Par­tons du prin­cipe que vous avez ins­tal­lé BLAST. Com­men­çons par pré­pa­rer la "que­ry" : les séquences CRISPR à recher­cher. Pour cela on crée un fichier "crispr.fasta", de ce type :

>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 :

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 :

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

[...]

>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 :

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 :
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 :
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 :
>>> 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 :
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 :
>>> 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 :

 

Merci à Sylvain P., Clem_, ook4mi et Yoann M. pour leur relecture.

 



Pour continuer la lecture :


Commentaires

9 réponses à “BLAST en pratique”

  1. Avatar de Jacques D
    Jacques D

    J’aime beau­coup cet article pour son esprit pra­tique. Le Blast et le “par­sage” de ses resul­tats sont des ele­ments recur­rent en bioin­fo. Il est donc bon de don­ner des fiches pra­tiques pour tout nou­veau uti­li­sa­teur (bio­in­for­ma­ti­cien et autres) ou rafrai­chir nos memoires trop sou­vent faillibles.

    Je sou­hai­te­rai contri­buer au sujet :
    — D’abord une infor­ma­tion utile pour les cas des courtes sequences comme celle de 20 nucleo­tides du cas pra­tique decrit ici. Le “Word size” par defaut, par exemple sur le site du ncbi, est de 28. Ce qui signi­fie qu’il sera impos­sible de trou­ver des sequences qui s’aligne sur moins de 28 nucleo­tides… Mais ils ont pense a tout : “Ins­tead, the nucleo­tide and pro­tein blast pro­grams auto­ma­ti­cal­ly check for short que­ries and adjust the search para­me­ters accor­din­gly. This adjust­ment occurs when the que­ry, either nucleo­tide or ami­no acid, is of length 30 or less.” Donc dans ce cas on aura bien un resul­tat pour cette sequence. Appa­re­ment la ver­sion de blast que tu uti­lise ici en local arbore le meme com­por­te­ment.
    La ver­sion “Nucleo­tide-Nucleo­tide 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 pen­ser a reduire le para­metre (max 15 … min 4) pour obte­nir un resul­tat.

    - Ensuite j’aime beau­coup uti­li­ser le bash pour se genre d’etude. Je vous pre­sente une maniere de pro­ce­der pour avoir les sequences pres des telo­meres qui sont absentes des regions inter­ge­no­miques :

    1) blastn ‑db reference.fa ‑que­ry crisp.fa ‑outfmt '6 qse­qid sse­qid pident length mis­match gapo­pen qstart qend sstart send slen' ‑word_​size 4 > resu.txt

    Cette com­mande ecrit les resul­tats du blast au for­mat tabule dans le fichier resu.txt. (Le for­mat tabule se choi­sit avec l’option 6. Vous pou­vez vous-meme choi­sir d’autres info ou chan­ger 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 telo­me­ric){ if (inter­ge­no­mic [i] == "") print i}}' resu.txt);
    for i in $seqs ; do
    awk ‑v test=$i '{if ($1 == test ) print $0}' resu.txt
    done

    La com­mande awk peut sem­bler com­plexe (pour rap­pel, $0=toute la ligne cou­rante, $1=colonne 1 de la ligne cou­rante, $2= colonne 2 de la ligne cou­rante, etc. Et bien sur il par­court les lignes une par une de haut en bas).
    Pour fil­trer les resul­tats de blasts awk est tres utile. Par exemple, dans mon resul­tat mon pour­cen­tage d’identite (pident) se retrouve dans la colonne 3. Pour affi­cher seule­ment les resul­tats de blast qui ont un pident > 95 il suf­fit de faire :
    awk '{if($3>95) print $0}' resu.txt
    on peut ain­si rajou­ter autant de test que neces­saire …

    P.S : Je sais pas si ca existe deja sur le site, mais il serait peut-etre bon d’ajouter aus­si l’approche en Bio­Perl…
    P.S2 : Desole j’ai pas d’accent sur mon cla­vier …

    1. Mer­ci pour ce com­men­taire !

      D'apres ce que j'ai lu, c'est Mega­blast qui a un word size de 28, et blastn est a 11 :

      "The most impor­tant rea­son that blastn is more sen­si­tive than MEGABLAST is that it uses a shor­ter default word size (11)". (http://​blast​.ncbi​.nlm​.nih​.gov/​B​l​a​s​t​.​c​g​i​?​C​M​D​=​W​e​b​&​P​A​G​E​_​T​Y​P​E​=​B​l​a​s​t​D​o​c​s​&​D​O​C​_​T​Y​P​E​=​P​r​o​g​S​e​l​e​c​t​i​o​n​G​u​i​d​e​#​t​a​b31, 4.2).

      Par contre c'est vrai qu'il faut bien pen­ser a l'ajuster car c'est le para­metre le plus impor­tant, par exemple un blastn avec des sequences de lon­gueur >100.

      Mer­ci aus­si pour la ver­sion bash. Je ne connais pas assez Perl mais je pour­rais ecrire une ver­sion R si quelqu'un s'y inter­esse.

      1. Avatar de Jacques D
        Jacques D

        C'est une bonne idée une ver­sion R !

        Autant pour moi, pour le "Word Size" tu as rai­son. Par défaut il te mette sur mega­blast quand tu choi­sis blastn sur NCBI … d'ou ma confu­sion.
        Par contre, il est inté­res­sant de noter que mega­blast a ce com­por­te­ment adap­ta­tif en fonc­tion de la lon­gueur de la séquence, dont je par­lais pré­cé­dem­ment (http://​blast​.ncbi​.nlm​.nih​.gov/​B​l​a​s​t​.​c​g​i​?​C​M​D​=​W​e​b​&​P​A​G​E​_​T​Y​P​E​=​B​l​a​s​t​D​o​c​s​&​D​O​C​_​T​Y​P​E​=​FAQ).
        Par contre je com­prends tou­jours pas pour­quoi ma ver­sion (tele­charge sur le NCBI) a un Word Size de 56 … le côte obs­cur de la bioin­fo sur­ement !

        P.S : les accents ont été ajoute par auto-com­ple­tion … pas ter­rible.

  2. Avatar de Mousset Carla
    Mousset Carla

    Je suis étu­diante en Mas­ter 2 en Eco­lo­gie et j'aimerai s'il vous plait savoir com­ment faire pour choi­sir la bonne séquence simi­laire résul­tat après avoir lan­cé une requête sur BLAST en ligne.

    Pour être plus pré­cise, j'ai des résul­tats de séquen­çage pour les­quels j'ai des séquences de nucléo­tide iden­ti­fiées à l'ordre par exemples, dont plu­sieurs pro­po­si­tions d'espèces ont été faites comme résul­tat de l'identification. Pour ces cas lorsque je blast la séquence j'ai effec­ti­ve­ment toutes ces pro­po­si­tions d'espèces tou­te­fois, je ne sais pas com­ment choi­sir la bonne séquence.

    Mer­ci d'avance pour votre réponse

    1. Bon­jour, mal­heu­reu­se­ment je ne me sou­viens plus du tout de tout ça, mais sur la ver­sion en ligne, je me sou­viens que les résul­tats appa­raissent triés de haut en bas du meilleur au moins bon, et ceci est pos­sible parce que chaque ali­gne­ment a un score, comme dans mon exemple où on peut lire " Score = 30.2 bits (15), Expect = 7.4".

  3. Avatar de coraline
    coraline

    Bon­jour,

    Je suis étu­diante en l2 génie bio­lo­gie option bio-infor­ma­tique, je me deman­dais s'il était pos­sible de lan­cer en local une com­mande per­met­tant d'obtenir la taille d'un frag­ment entre deux amorces.

    Je m'explique, je connais mes deux amorces, reward et reverse et j'ai télé­char­gé la séquence de la bac­té­rie qui m'interesse, je vou­drais donc pou­voir faire mat­cher mes deux amorces avec cette séquence_​ref et obte­nir la taille de la séquence entre les deux amorces (ici ça me per­met­trai d'otenir la taille de l'ARNr S16) dans un fichier result.
    Est-ce pos­sible et si oui pour­riez vous m'aider pour la com­mande ?

    Mer­ci d'avance pour votre réponse,
    Cor­dia­le­ment,
    Cora­line.

  4. Bon­jour ya t'il un moyen de bous contac­ter ? C'est au sujet de la géné­tique..

    1. Avatar de Yoann M.
      Yoann M.

      Bon­jour rk,

      Pour nous contac­ter nous vous invi­tons à nous envoyer un mail à l'adresse qui figure en bas de page d'un article : "Pour tout ren­sei­gne­ment sup­plé­men­taire : admin AT bioin​fo​-fr​.net"

  5. Avatar de Ismael

    Bon­jours
    Je suis étu­diant en mas­ter bio­chi­mie, je vou­drai savoir com­ment les résul­tats de blast sont-ils orga­ni­sés ?

Laisser un commentaire