- Le blog participatif de bioinformatique francophone depuis 2012 -

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 :

Ensuite, comme nous four­nis­sons nous-mêmes la réfé­rence "reference.fasta" (pour l'exemple, on pour­rait prendre le chro­mo­some 20 de l'humain : chr20.fa.gz, à décom­pres­ser et renom­mer), il faut créer l'index cor­res­pon­dant avec la com­mande "for­matdb" en ligne de com­mande :

Enfin on peut lan­cer BLAST en lui don­nant la que­ry et la réfé­rence, comme ça :

Le résul­tat, "alignments.txt", peut être assez lourd. Pour chaque séquence (par ex. chro­mo­some) de la réfé­rence, il donne tous les ali­gne­ments trou­vés de chaque séquence de crispr.fasta. Un extrait cor­res­pon­dant à la que­ry "seq3" (en uti­li­sant le chro­mo­some 20 de l'humain — hg19 — comme réfé­rence) :

>chr20
Length = 63025520

Score = 30.2 bits (15), Expect = 7.4
Iden­ti­ties = 15/​15 (100%)
Strand = Plus /​ Minus

Que­ry : 1 cgccccctcgtggcg 15
|||||||||||||||
Sbjct : 1946483 cgccccctcgtggcg 1946469

Score = 30.2 bits (15), Expect = 7.4
Iden­ti­ties = 15/​15 (100%)
Strand = Plus /​ Plus

Que­ry : 3 ccccctcgtggcgcc 17
|||||||||||||||
Sbjct : 34330497 ccccctcgtggcgcc 34330511

[…]
On voit que "seq3" a trou­vé deux homo­logues sur le chro­mo­some 20 : les bases 1–15 en posi­tion 1946483–1946469, et les bases 3–17 en 34330497–34330511. Un score d'enrichissement est don­né, ain­si que le nombre de matches sans erreur ("Iden­ti­ties"). Tout en haut est indi­quée la taille du chro­mo­some ; on voit ain­si que ces deux ali­gne­ments par­ti­cu­liers sont loin des télo­mères, donc pas du tout ce qui nous inté­resse.

Mais quand il y a des mil­liers de lignes comme celles-ci, on se dit qu'il faut un moyen plus sys­té­ma­tique d'examiner les résul­tats. De toute évi­dence, le par­ser soi-même est un cau­che­mar ; heu­reu­se­ment, des outils sont déjà dis­po­nibles. J'ai choi­si Bio­py­thon, la librai­rie bioin­fo du lan­gage Python. Voi­là com­ment exa­mi­ner le conte­nu du fichier :

A vous ensuite d'adapter le script à vos besoins par­ti­cu­liers.

Remarques :

Mer­ci à Syl­vain P., Clem_​, ook4mi et Yoann M. pour leur relec­ture.




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