BLAST

Prérequis

Cette section prolonge les sections 2.4 (alignement local de Smith-Waterman) et 2.5 (matrices de substitution PAM et BLOSUM) : vous devez savoir construire une matrice de substitution, attribuer une pénalité de gap et calculer un alignement local par programmation dynamique. Nous réutilisons aussi les formats de séquences de 1.2 et la lecture de fichiers avec Biopython de 1.3. Les calculs statistiques ne demandent que la bibliothèque standard de Python ; Biopython 1.88 sert uniquement à dialoguer avec NCBI.

Introduction

Comparer une séquence à une banque est l'opération la plus fréquente de la bioinformatique. Un généticien qui assemble un génome veut savoir à quelle espèce ressemble un contig ; un biologiste qui clone un gène veut identifier ses paralogues ; un épidémiologiste qui séquence un virus veut rattacher une souche à une lignée connue. Dans tous les cas, la question biologique est la même : existe-t-il dans la banque une séquence apparentée, et cette ressemblance est-elle plus forte que ce que le hasard produirait ? La réponse ne se lit pas dans un pourcentage d'identité isolé ; elle se construit avec un score, une longueur de recherche et une statistique.

Le défi computationnel est redoutable. L'algorithme local de Smith-Waterman trouve l'alignement optimal, mais son coût est quadratique en temps et en mémoire pour chaque paire. Comparer une protéine de 400 acides aminés à une banque de plusieurs centaines de millions de résidus demanderait des dizaines de milliards de cellules de programmation dynamique, répétées pour chaque requête. Même avec du matériel moderne, cette approche exhaustive est hors de portée pour une recherche interactive. BLAST (Basic Local Alignment Search Tool) accepte donc de renoncer à la garantie d'optimalité globale pour gagner plusieurs ordres de grandeur : il cherche d'abord des graines exactes courtes, puis étend les régions prometteuses.

Cette heuristique a deux conséquences. BLAST peut manquer des homologues divergents, et un score élevé ne prouve pas la parenté : il faut le confronter à la taille de l'espace de recherche. La E-value quantifie ce compromis. Nous construisons les briques du calcul, puis nous interrogeons NCBI avec Biopython.

Théorie

Pourquoi une heuristique plutôt qu'un alignement optimal

L'alignement local optimal maximise le score sur tous les couples de positions. Sa complexité $O(mn)$ vient de la table de programmation dynamique : chaque cellule combine le score diagonal et les extensions de gaps. Pour une requête de longueur $m$ et une banque de taille totale $n$, une recherche exhaustive coûte $O(mn)$ par requête, avec un facteur constant important. Une banque protéique de $2 \times 10^{8}$ résidus et une requête de 400 acides aminés donnent $8 \times 10^{10}$ cellules : impossible à recalculer pour chaque utilisateur.

BLAST inverse la stratégie : il repère des correspondances courtes, puis ne calcule que dans leurs voisinages. Le résultat n'est pas garanti optimal, mais il est généralement excellent pour des séquences apparentées, ce qui explique sa place comme outil de première intention.

Graines et mots de taille k

Une graine (seed) est une correspondance exacte entre un mot de longueur $k$ de la requête et un mot de la banque. Pour BLASTN, les mots sont des suites de nucléotides ; pour BLASTP, des suites d'acides aminés. On indexe les mots de la requête, puis on parcourt la banque : chaque mot partagé signale une diagonale candidate.

La taille $k$ gouverne le compromis sensibilité–vitesse. Un petit $k$ produit beaucoup de graines, donc beaucoup de fausses pistes, mais détecte des homologues divergents. Un grand $k$ réduit les graines fortuites et accélère la recherche, au risque de manquer une substitution qui a détruit le seul mot exact. Avec la requête ACGTAC et la cible ACGAAC, les 4-mères ne partagent rien, alors que le 3-mère ACG est commun.

Les valeurs par défaut reflètent ce compromis : BLASTN utilise généralement 11 nucléotides, BLASTP 3 acides aminés avec un seuil de score, car les substitutions conservatrices sont fréquentes. BLASTP accepte aussi les mots dont le score de substitution dépasse un seuil, et exige souvent deux mots compatibles sur la même diagonale avant d'étendre (règle des deux coups).

Extension sans gap, puis avec gaps

Chaque graine est étendue le long de la diagonale. On additionne les scores de substitution et on suit le meilleur score courant. L'extension s'arrête lorsque le score chute de plus d'un seuil $X$ sous le maximum : c'est le critère de coupure X (X-dropoff). Le segment obtenu est un HSP (high-scoring segment pair), une paire de segments locaux sans gap.

Les HSP prometteurs passent dans une seconde phase avec gaps : un alignement dynamique local restreint à une fenêtre autour de la graine, avec pénalités d'ouverture et d'extension. Cette phase coûte plus cher mais s'applique à peu de candidats. Un même couple requête–cible peut contenir plusieurs HSP, car des domaines différents ont parfois des histoires évolutives distinctes.

Filtres de complexité

Les régions de faible complexité — suites de A, répétitions d'un dinucléotide, segments riches en un acide aminé — produisent des graines par hasard et gonflent les scores. Sans filtre, une requête comme AAAAAAAAAAAA trouve des milliers de cibles sans intérêt. BLAST masque ces régions avant la recherche : SEG pour les protéines, DUST pour les nucléotides. Le masquage peut être doux (la région reste dans l'alignement mais n'est pas utilisée comme graine) ou dur (elle est remplacée par N ou X).

Le filtre améliore la spécificité mais peut supprimer un vrai domaine de faible complexité. Il faut donc savoir s'il est actif et le signaler. Un hit qui n'apparaît qu'après désactivation du filtre mérite un examen. La complexité n'est pas un critère de significativité : c'est un prétraitement.

E-value, bit score et espace de recherche

La théorie de Karlin-Altschul décrit la distribution des scores d'alignements locaux entre séquences aléatoires de même composition. Pour un score brut $S$, l'espérance du nombre de HSP fortuits d'au moins ce score est

où $m$ et $n$ sont les longueurs effectives, et $K$ et $\lambda$ dépendent du système de scores et de la composition. La longueur effective est inférieure à la longueur brute, car les extrémités offrent moins d'occasions de score. La E-value n'est pas une probabilité : c'est un nombre attendu de hits au moins aussi bons. Une valeur de 0,001 signifie qu'environ un millier de recherches produiraient un tel hit par hasard, pas « 0,1 % de chance que ce hit soit vrai ».

Le bit score normalise le score brut pour le rendre comparable entre systèmes de scores :

Avec cette définition, $E = mn\,2^{-S'}$. Le bit score ne dépend pas de la taille de la banque ; la E-value, si. Deux hits ne se comparent donc jamais par leur seule E-value s'ils viennent de banques de tailles différentes. Le bit score classe les alignements, la E-value décide de leur significativité dans un contexte donné.

Exemple travaillé

Prenons des paramètres pédagogiques explicites : $K = 0{,}1$, $\lambda = 0{,}3$, une requête effective $m = 100$, une banque effective $n = 1\,000$ et un score brut $S = 50$. La E-value vaut

Le bit score correspondant est

Avec la même requête dans une banque de $10^{6}$ résidus effectifs, la E-value devient $3{,}059$. Le score et le bit score n'ont pas changé, mais le hit n'est plus significatif au seuil de $10^{-5}$ : le nombre d'occasions de produire un score de 50 par hasard a été multiplié par mille. Si le score passe à 60, la E-value tombe à $1{,}52 \times 10^{-4}$ dans la petite banque et à $0{,}152$ dans la grande. La conclusion dépend du couple (score, espace de recherche).

Identité, couverture et significativité

Le pourcentage d'identité mesure la proportion de colonnes identiques. Il ne dit rien de la longueur de l'alignement ni du nombre de candidats testés. Un segment de 20 acides aminés identiques peut avoir une E-value médiocre dans une grande banque ; à l'inverse, un alignement de 300 résidus à 30 % d'identité peut être très significatif, car le signal cumulé est fort.

La couverture de requête, ou qcov, rapporte la longueur alignée de la requête à sa longueur totale. Un hit qui couvre 15 % de la requête peut correspondre à un domaine partagé, pas à un orthologue complet. Pour annoter, examinez E-value, bit score, identité, couverture et nombre de HSP. Le meilleur hit peut être un paralogue, une protéine voisine ou un contaminant.

Les programmes BLAST

BLASTX et TBLASTN traduisent la séquence nucléotidique dans les six cadres de lecture, trois par brin. Un hit absent en BLASTN peut apparaître en BLASTX si le gène contient des introns ou si le cadre est décalé. TBLASTX combine les deux traductions et coûte cher. PSI-BLAST, fondé sur des profils, étend la sensibilité aux homologues distants.

Banques, versions et reproductibilité

Les banques NCBI usuelles sont nt, nr, refseq_protein, swissprot et pdb. Leur contenu change en permanence : une E-value calculée contre nr en 2020 n'est pas comparable à une E-value calculée en 2026, car la taille et la composition ont changé.

Pour une analyse reproductible, enregistrez programme et version, nom et version de la banque, date de téléchargement, commande complète, paramètres de score et de filtrage, et fichier de sortie. Un tableau de hits sans version de banque n'est pas un résultat vérifiable.

Implémentation

Les calculs locaux ci-dessous sont exécutés avec Python 3.13.5 et Biopython 1.88 le 2026-11-01 ; les sorties sont vérifiées. Les appels BLAST réels ne sont pas exécutés dans cet environnement : ils sont signalés comme illustratifs.

Indexer les graines

kmer_seeds parcourt la séquence une fois et range chaque mot avec ses positions. La casse est normalisée en majuscules, les chevauchements sont conservés et les positions sont en base 0, conformément aux indices Python.

def kmer_seeds(seq, k):
    if k <= 0:
        raise ValueError("k must be a positive integer")
    seq = seq.upper()
    if k > len(seq):
        return {}
    seeds = {}
    for start in range(len(seq) - k + 1):
        word = seq[start:start + k]
        seeds.setdefault(word, []).append(start)
    return seeds
 
print(kmer_seeds("ACGTACGT", 3))

Sortie vérifiée :

{'ACG': [0, 4], 'CGT': [1, 5], 'GTA': [2], 'TAC': [3]}

Le coût est $O(n)$ en temps et, au pire, $O(n)$ en mémoire. Pour un alphabet de taille $\sigma$, il existe au plus $\sigma^{k}$ mots différents : la table reste petite pour $k=3$, mais devient creuse et coûteuse pour $k=11$ sur l'ADN.

Calculer une E-value et un bit score

La fonction evalue prend le score brut, la constante $K$ (nommée k dans le code), les longueurs effectives et $\lambda$. Le bit score se déduit de la même formule. Les paramètres ci-dessous sont pédagogiques : ils ne sont pas les valeurs de BLOSUM62.

import math
 
def evalue(score, k, m, n, lam):
    if k < 0:
        raise ValueError("K must be non-negative")
    if m <= 0 or n <= 0:
        raise ValueError("m and n must be positive")
    if lam <= 0:
        raise ValueError("lambda must be positive")
    return k * m * n * math.exp(-lam * score)
 
score = 50
value = evalue(score, 0.1, 100, 1000, 0.3)
bit = (0.3 * score - math.log(0.1)) / math.log(2)
print(round(value, 6))
print(round(bit, 3))
print(round(evalue(score, 0.1, 100, 1_000_000, 0.3), 6))

Sortie vérifiée :

0.003059
24.962
3.059023

Une E-value peut devenir si petite que math.exp renvoie 0.0 par sous-dépassement : cela signifie « plus petit que la précision du flottant », pas « significativité infinie ». Conservez le bit score et la couverture pour garder la force relative des alignements.

Lire une sortie tabulaire

La sortie outfmt 6 contient une ligne par HSP, avec des colonnes séparées par des tabulations. Les douze premières sont standard : pident est un pourcentage, length la longueur de l'alignement avec gaps, et les coordonnées sont 1-based et inclusives. Les colonnes supplémentaires vont dans extra.

BLAST_FIELDS = (
    "qseqid", "sseqid", "pident", "length", "mismatch", "gapopen",
    "qstart", "qend", "sstart", "send", "evalue", "bitscore",
)
INTEGER_FIELDS = ("length", "mismatch", "gapopen", "qstart", "qend", "sstart", "send")
FLOAT_FIELDS = ("pident", "evalue", "bitscore")
 
def parse_blast_tabular(text):
    hits = []
    for line_number, raw_line in enumerate(text.splitlines(), start=1):
        line = raw_line.strip()
        if not line or line.startswith("#"):
            continue
        fields = line.split("\t")
        if len(fields) < len(BLAST_FIELDS):
            raise ValueError(f"line {line_number}: expected at least 12 columns")
        hit = dict(zip(BLAST_FIELDS, fields[:12]))
        hit["extra"] = fields[12:]
        try:
            for name in INTEGER_FIELDS:
                hit[name] = int(hit[name])
            for name in FLOAT_FIELDS:
                hit[name] = float(hit[name])
        except ValueError as exc:
            raise ValueError(f"line {line_number}: invalid numeric field") from exc
        if not 0 <= hit["pident"] <= 100:
            raise ValueError(f"line {line_number}: pident outside [0, 100]")
        if hit["length"] <= 0:
            raise ValueError(f"line {line_number}: length must be positive")
        hits.append(hit)
    return hits

La commande BLAST+ correspondante, illustrative et non exécutée ici, serait :

blastp -query query.faa -db swissprot \
  -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" \
  -evalue 1e-5 -max_target_seqs 10 -out hits.tsv

Sur un tableau synthétique à identifiants fictifs, le parseur renvoie des types corrects :

sp|P00001|TEST1 99.5 200 1e-50 190.0
sp|P00002|TEST2 45.0 120 2e-10 60.5
sp|P00003|TEST3 28.0 80 0.02 35.2
float int float

Ces trois lignes sont une sortie illustrative : les identifiants TEST1 à TEST3 sont inventés pour l'exercice et ne correspondent à aucun hit réel. Elles servent uniquement à vérifier le découpage, les conversions et les filtres.

Interroger NCBI avec Biopython

Biopython encapsule l'interface Web de NCBI. Le code suivant, non exécuté ici car il dépend du réseau, montre la séquence typique : lire une requête FASTA, envoyer la recherche, conserver le XML, puis parser le résultat. La version testée de l'API est Biopython 1.88.

from Bio import SeqIO
from Bio.Blast import NCBIWWW, NCBIXML
 
record = SeqIO.read("query.faa", "fasta")
handle = NCBIWWW.qblast(
    "blastp",
    "swissprot",
    record.seq,
    expect=1e-5,
    hitlist_size=10,
    format_type="XML",
)
with open("blast_result.xml", "w", encoding="utf-8") as out:
    out.write(handle.read())
 
with open("blast_result.xml", encoding="utf-8") as out:
    for record in NCBIXML.parse(out):
        print(record.application, record.database)
        for alignment in record.alignments:
            for hsp in alignment.hsps:
                print(alignment.title, hsp.expect, hsp.bits, hsp.identities)

qblast attend la réponse de NCBI : une requête lente peut expirer. Prévoyez un délai, un nombre d'essais limité et un traitement des erreurs, et respectez les conditions d'utilisation de NCBI. format_type="XML" est indispensable pour NCBIXML ; une sortie tabulaire se traite avec le parseur précédent.

Dans le XML, record.database, record.database_sequences et record.database_length documentent la banque. Chaque alignment porte un titre, un identifiant et des HSP ; chaque HSP expose expect, bits, score, identities et align_length. Ces champs ne remplacent pas la lecture des coordonnées et de la couverture.

Étude de cas

Nous prenons la β-galactosidase LacZ d'Escherichia coli K-12 MG1655, portée par le chromosome NC_000913.3. Une recherche BLASTP contre une banque protéique retrouve-t-elle des β-galactosidases apparentées, et comment distinguer le signal principal du bruit de faible complexité ?

La démarche tient en cinq étapes. Extraire la séquence protéique annotée et noter l'accession du génome et la version de l'annotation. Lancer BLASTP 2.16 ou ultérieur contre Swiss-Prot, avec un seuil de $10^{-5}$, un nombre de cibles limité et SEG actif. Demander les colonnes standard, puis éventuellement qcovs. Conserver la sortie et la date de la banque. Enfin, filtrer et trier les hits avec un rapport qui sépare E-value, identité et couverture.

Les résultats attendus sont qualitatifs. Le meilleur hit doit être une β-galactosidase proche de LacZ, avec une E-value très inférieure au seuil, une identité élevée et une couverture proche de 100 %. Les suivants peuvent être des orthologues d'entérobactéries, puis des glycosidases plus distantes, avec une identité plus faible et une couverture partielle. Un hit court ou répétitif à E-value limite doit être examiné séparément. Aucun chiffre précis n'est avancé : il dépendrait de la version de Swiss-Prot.

Pour illustrer le traitement, voici un tableau synthétique à identifiants fictifs, choisis pour montrer le filtrage ; ce n'est pas un résultat BLAST réel.

# sortie illustrative, identifiants fictifs
q1	sp|P00001|TEST1	99.500	200	1	0	1	200	5	204	1e-50	190
q1	sp|P00002|TEST2	45.000	120	60	4	10	129	2	121	2e-10	60.5
q1	sp|P00003|TEST3	28.000	80	55	6	1	80	1	80	0.02	35.2

Avec une longueur de requête de 200, le premier hit couvre 100 % de la requête, le deuxième 60 %, le troisième 40 %. Un rapport exigeant une E-value inférieure ou égale à $10^{-5}$ et une identité d'au moins 30 % retient les deux premiers. Le meilleur hit a une E-value de $10^{-50}$ et un bit score de 190 ; la couverture confirme un alignement quasi complet.

Un hit à 99,5 % sur 200 acides aminés est probablement la même protéine ou un orthologue proche ; un hit à 45 % sur 120 acides aminés peut être un paralogue. La E-value ne tranche pas entre ces hypothèses : elle dit seulement que ces alignements sont improbables par hasard. La conclusion doit citer la banque, sa version et la date.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-2/2.6-blast -q

python -m pytest exercises/module-2/2.6-blast -q

python -m pytest exercises/module-2/2.6-blast/projet -q

Résumé

• BLAST est une heuristique seed-and-extend : des mots exacts ou proches amorcent des extensions locales, d'abord sans gap, puis avec gaps.
• La taille des mots $k$ arbitre entre sensibilité et vitesse ; un petit $k$ trouve plus d'homologues divergents mais produit plus de faux candidats.
• Les filtres SEG et DUST masquent les régions de faible complexité, dont les scores sont souvent fortuits ; leur état doit être documenté.
• La E-value $E = Kmn e^{-\lambda S}$ est le nombre attendu de hits fortuits au moins aussi bons, pas une probabilité de parenté.
• Le bit score $S' = (\lambda S - \ln K)/\ln 2$ ne dépend pas de la taille de la banque ; la E-value, si.
• Identité et couverture sont des observations ; la significativité se juge avec la E-value, le bit score et l'espace de recherche.
• BLASTN, BLASTP, BLASTX, TBLASTN et TBLASTX se choisissent selon la nature des séquences et la présence de cadres de lecture.
• Biopython 1.88 permet d'appeler qblast et de parser NCBIXML, à condition de conserver la sortie, la version de la banque et la date.

Références

• Altschul, S. F., Gish, W., Miller, W., Myers, E. W. & Lipman, D. J. (1990). Basic local alignment search tool. Journal of Molecular Biology, 215(3), 403–410.
• Karlin, S. & Altschul, S. F. (1990). Methods for assessing the statistical significance of molecular sequence features by using general scoring schemes. Proceedings of the National Academy of Sciences, 87(6), 2264–2268.
• Altschul, S. F. et al. (1997). Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Research, 25(17), 3389–3402.
• Camacho, C. et al. (2009). BLAST+: architecture and applications. BMC Bioinformatics, 10, 421.
• Cock, P. J. A. et al. (2009). Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics, 25(11), 1422–1423.
• NCBI. BLAST Documentation. https://blast.ncbi.nlm.nih.gov/doc/blast-help/
• NCBI. BLAST Command Line Applications User Manual. https://www.ncbi.nlm.nih.gov/books/NBK279690/