Contrôle qualité NGS

Prérequis

Cette section suppose acquis le format FASTQ et le décodage des scores Phred de 1.2 : quatre lignes par enregistrement, qualité de même longueur que la séquence, encodage Phred+33. Vous devez savoir calculer $P = 10^{-Q/10}$ et reconnaître qu'une moyenne de qualité n'est pas une probabilité. Le code n'utilise que la bibliothèque standard de Python 3.11+ ; FastQC et fastp sont présentés par leurs commandes et leurs rapports, sans être requis pour les exercices.

Introduction

Le séquençage de nouvelle génération (next-generation sequencing, NGS) produit des millions de lectures courtes, chacune accompagnée d'un score de confiance par base. Ces lectures ne sont pas des observations parfaites : la qualité chute souvent en fin de lecture, les adaptateurs de préparation de banque peuvent être lus lorsque l'insert est court, certaines régions sont sous-représentées pour des raisons de composition en GC, et des fragments peuvent être dupliqués par amplification. Un contrôle qualité mal conduit laisse passer des artefacts qui deviennent, plus loin, des variants faux, des couvertures déséquilibrées ou des pics d'expression trompeurs.

L'enjeu computationnel est double. Il faut mesurer, base par base et lecture par lecture, des propriétés qui ne tiennent pas dans un seul chiffre : une fraction Q30 globale peut être excellente alors qu'une minorité de lectures est inutilisable. Il faut ensuite décider quoi corriger, car le trimming retire des bases et le filtering des lectures entières ; les deux modifient les résultats en aval. Un pipeline reproductible conserve les paramètres, les versions et les rapports, pas seulement les fichiers nettoyés.

Nous suivons le fil rouge de Escherichia coli K-12 MG1655, dont le génome de référence est NC_000913.3, mais les principes valent pour tout organisme. Le jeu de démonstration est un petit FASTQ synthétique explicitement étiqueté : il ne prouve rien sur une banque réelle, il rend les calculs vérifiables à la main.

Théorie

Ce que contrôle réellement le contrôle qualité

Le contrôle qualité répond à une question limitée : ce jeu de lectures permet-il l'analyse prévue, dans des conditions documentées ? Il ne répond pas à la question biologique. Une lecture peut être de haute qualité et provenir d'une contamination ; une lecture de qualité moyenne peut contenir un vrai variant. La qualité Phred estime la probabilité qu'une base appelée soit fausse, pas la probabilité qu'elle soit biologiquement pertinente.

Les défauts se répartissent en trois familles : séquençage (qualité par base, N, poly-G, chute en 3'), banque (adaptateurs, dimères, composition GC, taille des inserts, duplication PCR) et échantillon (contamination, ARN ribosomal, ADN mitochondrial). FastQC observe des symptômes ; il ne nomme pas toujours la cause.

Qualité par base, Q30 et erreurs attendues

Le score Phred $Q$ d'une base vérifie $Q = -10\log_{10}P$, où $P$ est la probabilité d'erreur estimée. Q20 correspond à $P = 10^{-2}$, Q30 à $P = 10^{-3}$. La fraction Q30 d'un jeu de $N$ bases est

où la fonction indicatrice vaut 1 si la condition est vraie et 0 sinon. Cette fraction est pondérée par les bases, pas par les lectures : une lecture de 250 bases pèse plus qu'une lecture de 50 bases.

Le nombre attendu d'erreurs d'une lecture est $\sum_j P_j$, et la probabilité qu'au moins une base soit fausse vaut $1-\prod_j(1-P_j)$ sous indépendance. Une lecture de 150 bases toutes à Q30 a encore une probabilité d'au moins une erreur de $1-(1-0{,}001)^{150} \approx 0{,}139$. C'est pourquoi un Q30 élevé ne signifie pas « lecture parfaite ».

Prenons une lecture de 12 bases dont 8 à Q40 et 4 à Q2. La fraction Q30 vaut $8/12 \approx 0{,}667$. Le nombre attendu d'erreurs est $8 \times 10^{-4} + 4 \times 10^{-0{,}2} \approx 2{,}525$. La probabilité d'au moins une erreur vaut $1-(1-10^{-4})^8(1-0{,}630957)^4 \approx 0{,}981$. La lecture est donc presque certainement erronée au moins une fois, malgré huit bases excellentes. Une moyenne de $Q$ à 27,3 ne le montrerait pas clairement.

Longueur des lectures et qualité en 3'

Un protocole 2 × 150 pb doit produire des lectures d'environ 150 pb, avec une dispersion faible. Une queue de lectures courtes peut signaler un trimming antérieur, une dégradation de l'ARN, un insert court ou un problème de cycles. La longueur brute inclut les adaptateurs et les bases de faible qualité ; la longueur après trimming ne mesure donc pas nécessairement l'insert biologique.

La qualité diminue fréquemment en 3' à cause du phasage et du pré-phasage des cycles Illumina. Un trimming des queues de faible qualité est souvent justifié, mais il faut éviter de couper toute la lecture : une lecture trop courte ne s'aligne plus de façon unique, et BWA-MEM exige une longueur minimale de graine. Couper 20 bases partout peut aussi supprimer un signal biologique en fin de transcrit.

Adaptateurs, read-through et poly-G

L'adaptateur est une séquence synthétique ligaturée aux extrémités du fragment lors de la préparation de banque. Si l'insert est plus court que la longueur lue, le séquenceur lit l'adaptateur et éventuellement le début de l'insert suivant : c'est le read-through. La séquence de référence d'Illumina TruSeq commence souvent par AGATCGGAAGAGC. Un adaptateur peut aussi apparaître seul, sous forme de dimère, lorsque la ligature a échoué.

La détection repose sur trois signaux : une séquence surreprésentée en 3', une augmentation de la teneur en adaptateur le long de la lecture, et une longueur d'insert anormalement courte. Un simple find de l'adaptateur complet rate les adaptateurs partiels. Un chevauchement de trois bases apparaît par hasard dans une fraction non négligeable des lectures ; un seuil de quatre à six bases est un compromis courant.

Le poly-G est un artefact des plateformes Illumina à deux couleurs, notamment NextSeq et NovaSeq : lorsqu'aucun signal n'est détecté, l'appel de base produit des G successifs. fastp peut retirer les queues poly-G ; il faut distinguer cet artefact d'une vraie région riche en G, rare mais possible. Le poly-G n'est pas un adaptateur et ne se traite pas avec la même option.

Biais GC, duplication et séquences surreprésentées

La fraction GC d'une lecture est $f_{GC} = (n_G+n_C)/L$. La distribution des $f_{GC}$ doit refléter la composition de l'échantillon. Un décalage vers les faibles GC peut indiquer une contamination AT-riche ou un biais d'amplification ; un décalage vers les forts GC, une contamination GC-riche ou une sélection de fragments. Un profil bimodal évoque un mélange d'espèces, mais une banque métagénomique peut être multimodale.

La duplication mesure la répétition de fragments. Si $N$ fragments sont lus et $U$ sont uniques, le taux de duplication vaut $1-U/N$. Une partie est biologique, une partie vient de la PCR ; le taux attendu dépend du protocole (génome entier, RNA-seq, amplicon). Supprimer les duplicats sans réfléchir peut donc détruire le signal. Pour l'appel de variants, on préfère souvent les marquer plutôt que les supprimer.

Une séquence surreprésentée est une séquence identique dont la fréquence $f_s = c_s/N$ dépasse un seuil. Elle peut être un adaptateur, un dimère, un ARN ribosomal, une séquence PhiX utilisée comme contrôle, un primer ou un contaminant. FastQC signale ces séquences sans les identifier ; une recherche par similarité ou un classifieur taxonomique est nécessaire pour conclure. Un petit jeu de lectures produit facilement des surreprésentations par hasard : deux lectures identiques sur quatre représentent déjà 50 %.

Lectures appariées

En séquençage apparié, les deux lectures d'une paire proviennent des deux extrémités du même fragment. Les identifiants sont reliés (/1 et /2, ou un champ 1:N:0:), et les deux fichiers doivent rester dans le même ordre. Les qualités des deux mates sont corrélées mais pas identiques. Filtrer un seul mate casse la paire ; la plupart des outils retirent la paire entière ou écrivent les lectures orphelines dans des fichiers séparés.

L'appariement aide aussi la correction. Lorsque les deux lectures se chevauchent, un désaccord peut être corrigé en utilisant le mate de meilleure qualité. fastp active cette correction par défaut pour les données appariées et permet de la désactiver. Elle modifie les bases : il faut la documenter, surtout pour l'étude de variants à faible fraction allélique.

FastQC : les modules et leur lecture

FastQC 0.12.1 regroupe onze modules : Basic Statistics (nombre, longueur, GC, encodage), Per base sequence quality (qualité par position), Per tile sequence quality (défauts de flow cell), Per sequence quality scores (distribution des moyennes), Per base sequence content (A, C, G, T par position), Per sequence GC content, Per base N content, Sequence Length Distribution, Sequence Duplication Levels, Overrepresented sequences et Adapter Content.

Les statuts vert, jaune et rouge sont des heuristiques calculées module par module. Un rouge n'est pas un échec automatique ; un vert ne garantit rien, et un petit fichier produit des avertissements instables. Lisez Basic Statistics, la qualité par base, les adaptateurs et les séquences surreprésentées, puis la duplication et le GC selon le protocole.

Un signal classique combine une chute de qualité en 3' et une remontée de l'adaptateur. Un biais de composition dans les premières bases vient souvent de l'amorçage aléatoire ; un pic de GC décalé peut être réel ou révéler un mélange. Aucun de ces signaux ne donne directement une commande de trimming.

fastp : corrections en une passe

fastp 0.23.4 lit les lectures, retire les adaptateurs, filtre la qualité, coupe les extrémités, corrige les chevauchements et produit un rapport HTML et JSON. Pour les données appariées, la détection d'adaptateur par chevauchement est activée avec --detect_adapter_for_pe ; pour les données simple brin, il faut --detect_adapter_for_se ou fournir --adapter_sequence. Les contrôles de qualité utilisent --qualified_quality_phred et --unqualified_percent_limit : une lecture est écartée si trop de ses bases sont sous le seuil. La longueur minimale se règle avec --length_required.

Le trimming des extrémités se fait avec --cut_front, --cut_tail ou --cut_right ; le poly-G avec --trim_poly_g, souvent activé par défaut sur les plateformes détectées, et désactivable avec --disable_trim_poly_g. La correction de chevauchement s'active par défaut en apparié et se désactive avec --disable_correction. Chaque option change le jeu de lectures : conservez la commande exacte, la version et les rapports avec les données.

Trimming et filtering : deux opérations différentes

Le trimming enlève des bases : adaptateur, queue de faible qualité, bases poly-G. La lecture reste présente, plus courte. Le filtering enlève des lectures entières : trop courtes, trop de bases sous le seuil, trop de N, mate non apparié. fastp fait les deux dans une seule passe.

Le trimming préserve la profondeur mais peut biaiser la longueur et supprimer du signal en 3'. Le filtering préserve la longueur des lectures restantes mais réduit la profondeur et peut créer un biais de sélection si les lectures écartées ne sont pas aléatoires. Un trimming agressif peut faire disparaître la couverture des extrémités de gènes ou les variants proches des extrémités. Fixez les seuils avant de regarder les résultats, puis vérifiez l'effet sur la couverture et l'appel de variants.

MultiQC : agréger sans masquer

MultiQC lit les rapports de FastQC, fastp, samtools, BWA et d'autres outils, puis produit un rapport HTML unique par échantillon et par métrique. Il agrège des sorties existantes sans refaire les calculs. Vérifiez que tous les échantillons sont présents : un rapport manquant est parfois silencieux.

MultiQC est utile pour comparer des dizaines de banques. Il ne remplace pas la lecture d'un rapport individuel : une moyenne de Q30 sur dix échantillons peut cacher un échantillon catastrophique.

Limites : ne pas confondre qualité et signal biologique

La qualité Phred est une estimation du base caller, pas une vérité biologique. Une base à Q40 peut être fausse, une base à Q20 peut être vraie. Une séquence surreprésentée peut être un transcrit abondant, un ARN ribosomal ou un primer.

Le trimming excessif est un risque réel. Couper les lectures à Q30 partout peut supprimer des régions réellement couvertes, réduire la diversité des inserts et fausser l'estimation de la taille des fragments. À l'inverse, ne rien couper peut laisser des adaptateurs qui seront soft-clippés par l'aligneur, ou des queues erronées qui deviendront des variants. La bonne décision dépend de l'objectif : appel de variants germinaux, détection somatique à faible fraction, quantification RNA-seq ou assemblage.

Enfin, le contrôle qualité ne détecte pas tout. Une contamination croisée à faible taux, un échange d'échantillons ou un biais de lot peuvent passer inaperçus. Les métadonnées, les témoins et la comparaison entre échantillons restent indispensables.

Implémentation

Les quatre fonctions ci-dessous sont volontairement courtes et sans dépendance. Elles illustrent les calculs que FastQC et fastp réalisent à grande échelle : décodage Phred+33, fraction Q30, trimming d'adaptateur, distribution des longueurs et détection de séquences identiques. Le corrigé complet est dans exercises/module-3/3.1-controle-qualite-ngs/solution_3_1.py.

PHRED33_OFFSET = 33
MAX_PHRED33 = 93
MIN_PARTIAL_OVERLAP = 4
 
 
def q30_fraction(qualities):
    total = 0
    passing = 0
    for quality in qualities:
        if not isinstance(quality, str):
            raise TypeError("each quality must be a string")
        for char in quality:
            score = ord(char) - PHRED33_OFFSET
            if score < 0 or score > MAX_PHRED33:
                raise ValueError(f"invalid Phred+33 quality character {char!r}")
            total += 1
            if score >= 30:
                passing += 1
    return passing / total if total else 0.0
 
 
def trim_adapters(seq, adapter):
    if not adapter:
        raise ValueError("adapter must not be empty")
    upper_seq = seq.upper()
    upper_adapter = adapter.upper()
    position = upper_seq.find(upper_adapter)
    if position != -1:
        return seq[:position]
    max_overlap = min(len(upper_adapter), len(upper_seq))
    for overlap in range(max_overlap, MIN_PARTIAL_OVERLAP - 1, -1):
        if upper_seq.endswith(upper_adapter[:overlap]):
            return seq[:-overlap]
    return seq
 
 
def read_length_distribution(reads):
    distribution = {}
    for read in reads:
        sequence = read["sequence"] if isinstance(read, dict) else read
        distribution[len(sequence)] = distribution.get(len(sequence), 0) + 1
    return dict(sorted(distribution.items()))
 
 
def detect_overrepresented(reads, min_fraction):
    if not 0.0 <= min_fraction <= 1.0:
        raise ValueError("min_fraction must be between 0 and 1")
    sequences = [
        (read["sequence"] if isinstance(read, dict) else read).upper()
        for read in reads
    ]
    total = len(sequences)
    if total == 0:
        return []
    counts = {}
    for sequence in sequences:
        counts[sequence] = counts.get(sequence, 0) + 1
    hits = [
        {"sequence": sequence, "count": count, "fraction": count / total}
        for sequence, count in counts.items()
        if count / total >= min_fraction
    ]
    hits.sort(key=lambda hit: (-hit["count"], hit["sequence"]))
    return hits

Les appels suivants donnent des résultats déterministes, vérifiés avec Python 3.13 :

ADAPTER = "AGATCGGAAGAGC"
qualities = ["IIIIIIIIIIII", "IIIIIIII####"]
print(round(q30_fraction(qualities), 4))
print(trim_adapters("ACGTACGTACGT" + ADAPTER, ADAPTER))
print(trim_adapters("ACGTACGTAGATCG", ADAPTER))
print(read_length_distribution(["ACGT", "ACGTACGT", "TTTT"]))
print(detect_overrepresented(["ACGT", "ACGT", "TTTT", "GGGG"], 0.5))

0.8333
ACGTACGTACGT
ACGTACGT
{4: 2, 8: 1}
[{'sequence': 'ACGT', 'count': 2, 'fraction': 0.5}]

Deux points sont critiques : q30_fraction valide les caractères, et trim_adapters accepte un chevauchement partiel d'au moins quatre bases pour les lectures qui s'arrêtent au milieu de l'adaptateur. Ce seuil limite les faux positifs sans les éliminer.

Pour un vrai jeu de lectures, FastQC et fastp s'exécutent en ligne de commande. Les commandes ci-dessous sont données avec les versions du cours ; leurs sorties ne sont pas reproduites ici et seraient illustratives.

fastqc --version
fastqc -t 4 -o qc/ reads_R1.fastq reads_R2.fastq
 
fastp \
  -i reads_R1.fastq -I reads_R2.fastq \
  -o reads_R1.trim.fastq -O reads_R2.trim.fastq \
  --detect_adapter_for_pe \
  --qualified_quality_phred 20 \
  --unqualified_percent_limit 20 \
  --length_required 50 \
  --trim_poly_g \
  --correction \
  --html fastp.html --json fastp.json --thread 4
 
multiqc qc/ fastp.html -o multiqc

fastp écrit deux fichiers synchronisés : les mates restent dans le même ordre. --length_required 50 filtre les lectures trop courtes après trimming, et les deux options de qualité écartent une lecture si plus de 20 % de ses bases sont sous Q20. Ces valeurs sont un point de départ.

Étude de cas

Nous travaillons sur le fichier fourni datasets/ngs/reads_R1.fastq, un exemple synthétique de deux lectures appariées avec reads_R2.fastq. Chaque lecture mesure 12 pb. La seconde lecture contient une queue de quatre bases à Q2, ce qui permet d'observer une chute de qualité localisée.

La question est : ce jeu est-il propre, et quelles corrections seraient justifiées ? Nous calculons la fraction Q30 sur les 24 bases de R1. Les 12 bases de pair1/1 et les 8 premières de pair2/1 sont à Q40 ; les 4 dernières de pair2/1 sont à Q2. La fraction Q30 vaut donc $20/24 \approx 0{,}833$. Le nombre attendu d'erreurs est $20 \times 10^{-4} + 4 \times 10^{-0{,}2} \approx 2{,}526$. Pour pair2/1, la probabilité d'au moins une erreur est $1-(1-10^{-4})^8(1-0{,}630957)^4 \approx 0{,}981$.

Le verdict global est trompeur : 83,3 % de bases Q30 semble acceptable, mais une lecture sur deux est presque certainement erronée. La moyenne de qualité vaut 33,7, ce qui masque encore la queue à Q2. Il faut lire la qualité par position. Un trimming de la queue de pair2/1 est raisonnable si le protocole exige Q20 ; filtrer la paire entière est plus agressif, mais cohérent si aucune base à Q2 n'est tolérée.

La détection d'adaptateur sur ces deux lectures renvoie zéro : la séquence AGATCGGAAGAGC n'apparaît ni dans pair1/1 ni dans pair2/1. Deux lectures ne suffisent pas à estimer un taux faible. Pour illustrer le read-through, ajoutons un enregistrement synthétique de 25 pb, ACGTACGTACGT suivi de l'adaptateur complet, avec 25 bases à Q40. trim_adapters trouve l'adaptateur à la position 0-based 12 et renvoie les 12 pb d'insert. Sans trimming, l'aligneur devrait gérer 13 pb d'adaptateur, soit un soft-clip ou un alignement trompeur.

read = "ACGTACGTACGT" + "AGATCGGAAGAGC"
trimmed = trim_adapters(read, "AGATCGGAAGAGC")
print(len(read), len(trimmed), trimmed)

25 12 ACGTACGTACGT

L'interprétation doit rester prudente : le fichier est synthétique et minuscule. Sur un jeu réel, on confronterait qualité par base, adaptateurs, inserts et duplication. On vérifierait aussi que le mate R2 présente la même chute : c'est le cas ici, ce qui oriente vers un problème de cycles plutôt que vers un défaut propre à une lecture.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-3/3.1-controle-qualite-ngs -q

python -m pytest exercises/module-3/3.1-controle-qualite-ngs -q

python -m pytest exercises/module-3/3.1-controle-qualite-ngs/projet -q

Résumé

• Le contrôle qualité répond à une question d'usage : le jeu de lectures permet-il l'analyse prévue, avec des paramètres documentés ?
• La fraction Q30 est pondérée par les bases ; une moyenne de qualité ou un Q30 global peut cacher une queue catastrophique et une lecture presque certainement fausse.
• Les adaptateurs produisent un read-through lorsque les inserts sont courts ; détectez aussi les chevauchements partiels et traitez les deux mates ensemble.
• fastp combine trimming, filtering, correction de chevauchement et détection d'adaptateurs ; conservez les paramètres et les rapports.
• FastQC décrit des symptômes par module ; ses statuts sont des heuristiques, pas des verdicts. MultiQC agrège des rapports existants sans remplacer leur lecture.
• Un mauvais signal de qualité n'est pas forcément un défaut biologique : duplication, GC extrême et surreprésentation peuvent être attendus selon le protocole.

Références

• Andrews, S. (2010). FastQC: A Quality Control Tool for High Throughput Sequence Data. Babraham Bioinformatics. https://www.bioinformatics.babraham.ac.uk/projects/fastqc/
• Chen, S., Zhou, Y., Chen, Y. & Gu, J. (2018). fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics, 34(17), i884–i890.
• Ewels, P., Magnusson, M., Lundin, S. & Käller, M. (2016). MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics, 32(19), 3047–3048.
• Babraham Bioinformatics. FastQC documentation (version 0.12.1). https://www.bioinformatics.babraham.ac.uk/projects/fastqc/
• OpenGene. fastp documentation (version 0.23.4). https://github.com/OpenGene/fastp
• MultiQC. MultiQC documentation. https://multiqc.info/docs/