Formats FASTA et FASTQ

Prérequis

Cette section prolonge 1.1 : vous devez savoir qu'une séquence est une chaîne sur un alphabet fini, que sa longueur se compte après concaténation des lignes et que les coordonnées biologiques sont 1-based. Les scripts n'utilisent que la bibliothèque standard de Python 3.11+ : str, gzip, math. Aucune connaissance de Biopython n'est nécessaire ; la bibliothèque sera présentée plus loin comme référence.

Introduction

Le séquençage produit des lectures (reads), c'est-à-dire des fragments courts lus sur un brin d'ADN. Pour exploiter ces lectures, il faut les stocker avec deux informations indissociables : la suite de nucléotides et la confiance que l'on accorde à chaque base. Le format FASTA ne conserve que la séquence ; le format FASTQ ajoute une qualité par base. On les trouve dans les génomes de référence, les banques de protéines et les fichiers de lectures brutes. Une erreur de lecture ou d'interprétation se propage ensuite dans l'alignement, l'appel de variants, le comptage d'expression et l'annotation.

Le défi computationnel est trompeur. Un fichier FASTA ressemble à quelques lignes de texte, et un fichier FASTQ à quatre lignes répétées. Pourtant, les fichiers réels contiennent des en-têtes sur plusieurs lignes, des lignes vides, des fins de ligne CRLF, des identifiants dupliqués, des qualités de longueur incohérente et parfois des archives gzip. Un parseur qui fonctionne sur un fichier propre peut échouer silencieusement sur un fichier téléchargé ailleurs. La qualité est particulièrement fragile : une chaîne ASCII décalée d'un caractère ne change pas la longueur du fichier, mais change toutes les probabilités d'erreur.

Nous partons du fil rouge du cours, Escherichia coli K-12 MG1655 (NC_000913.3), pour un extrait FASTA réel. Nous construisons aussi un petit FASTQ synthétique, étiqueté comme tel, qui contient une lecture de mauvaise qualité. L'objectif est de lire, vérifier, écrire et interpréter ces fichiers avec un code court, explicite et testable.

Théorie

La structure d'un fichier FASTA

Un fichier FASTA contient un ou plusieurs enregistrements. Chaque enregistrement commence par un en-tête > en première colonne ; le premier mot est l'identifiant, le reste une description libre. Les lignes suivantes contiennent la séquence, éventuellement coupées en lignes de 60 ou 80 caractères. Formellement, un enregistrement est un couple $(h, s)$ avec $s \in \Sigma^{*}$ ; sa longueur est $|s|$, pas le nombre de lignes.

>NC_000913.3:1-60 Escherichia coli str. K-12 substr. MG1655, complete genome
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTC

L'identifiant est NC_000913.3:1-60, la description commence à Escherichia. NCBI recommande un SeqID unique, sans espace, et des lignes d'au plus 80 caractères. Un enregistrement peut avoir une séquence vide, des identifiants dupliqués et des minuscules, souvent signe de soft-masking.

Les espaces et tabulations internes sont retirés par la plupart des parseurs : utile pour les alignements, mais cela peut masquer une corruption. Un FASTA peut contenir des symboles IUPAC ambiguus (N, R, Y) et, pour les protéines, * pour un codon stop. La validation de l'alphabet est une étape distincte de l'analyse syntaxique.

Le format FASTA original tolérait des lignes de commentaire commençant par ;, mais cette convention n'est pas universelle : Biopython ignore le texte avant le premier enregistrement et traite une ligne ; interne comme de la séquence. Notre parseur ignore ; avant le premier en-tête et lève une erreur à l'intérieur d'un enregistrement : mieux vaut refuser un fichier ambigu que produire une séquence fausse.

La structure d'un fichier FASTQ

Un enregistrement FASTQ moderne tient sur quatre lignes : une ligne d'en-tête qui commence par @, la séquence, une ligne qui commence par +, puis une chaîne de qualité ASCII de même longueur que la séquence. Formellement, un enregistrement est un triplet $(h, s, q)$ avec $|q| = |s|$. La ligne + peut répéter l'en-tête, mais cette répétition est facultative et souvent ignorée.

@read1
ACGTACGTACGT
+
IIIIIIIIIIII

La qualité n'est pas un nombre séparé : c'est une chaîne de caractères, où chaque caractère encode un score Phred. Le caractère I a le code ASCII 73 ; avec l'encodage Phred+33, il représente $Q = 73 - 33 = 40$. Le caractère # a le code 35, donc $Q = 2$. Une qualité peut contenir @ et + : c'est pourquoi il ne faut jamais découper un FASTQ sur ces caractères. La structure est positionnelle, pas lexicale.

La variante à quatre lignes est celle des fichiers Illumina modernes. Des fichiers anciens peuvent enrouler la séquence et la qualité sur plusieurs lignes ; un parseur complet doit alors suivre un état. Notre parseur traite la variante moderne et signale une longueur incohérente.

Les lignes vides entre deux enregistrements sont tolérées par beaucoup d'outils. À l'intérieur d'un enregistrement, une ligne vide est ambiguë : lecture vide ou ligne supprimée par un éditeur. Il faut la signaler. Enfin, un FASTQ peut être compressé en gzip ; le fichier doit être lu en mode texte après décompression.

Scores Phred, probabilités et encodages

Le score Phred $Q$ d'une base est défini par

où $P$ est la probabilité estimée que cette base soit fausse. La réciproque est

Q20 correspond à $P = 0{,}01$ (une erreur sur 100), Q30 à $P = 0{,}001$ (une sur 1 000). Les seuils Q20/Q30 comptent la proportion de bases au-dessus du seuil. Un Q30 n'est pas une garantie biologique : c'est une estimation du base caller.

L'encodage ASCII ajoute un décalage. En Phred+33 (Sanger, Illumina depuis CASAVA 1.8), ! vaut Q0 et ~ vaut Q93. En Phred+64 (Illumina 1.3–1.7), I vaut Q9, pas Q40. Les fichiers Solexa utilisent un autre système, fondé sur des log-odds. Un fichier dont tous les caractères sont au-dessus de @ peut relever des deux encodages : la détection est ambiguë, consultez les métadonnées.

Pour une lecture, la probabilité qu'au moins une base soit fausse n'est pas la moyenne des $Q$. En supposant les erreurs indépendantes,

et le nombre attendu d'erreurs est $\sum_{j=1}^{n} P_j$. La moyenne des scores, $\bar{Q}$, est une statistique descriptive commode, mais elle ne se convertit pas directement en probabilité.

Exemple travaillé à la main

Prenons la chaîne de qualité IIII#, de longueur 5. Les codes ASCII sont 73, 73, 73, 73 et 35. En Phred+33, les scores sont donc $Q = 40, 40, 40, 40, 2$. Les probabilités d'erreur par base valent $10^{-4}$ pour chaque I et $10^{-0.2} \approx 0{,}630957$ pour le #.

La probabilité que les cinq bases soient toutes correctes, sous l'hypothèse d'indépendance, est

La probabilité qu'au moins une base soit fausse est donc $1-0{,}368895 \approx 0{,}631105$, soit environ 63 %. Le nombre attendu d'erreurs vaut $4 \times 10^{-4} + 0{,}630957 \approx 0{,}631357$. Pourtant, la moyenne des scores est $(40+40+40+40+2)/5 = 32{,}4$, un chiffre qui donnerait une impression de très bonne qualité. C'est le piège central de cette section : une seule base catastrophique suffit à rendre la lecture peu fiable, et la moyenne ne le montre pas.

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$ : chaque base a 0,1 % de chance d'être fausse.

Complexité, mémoire et streaming

Pour $L$ caractères et $R$ enregistrements, un parseur coûte $O(L)$ en temps. Matérialiser tous les enregistrements coûte $O(L)$ en mémoire. En streaming, la mémoire supplémentaire est $O(\max_i |s_i|)$ pour la lecture courante ; les compteurs Q20/Q30 tiennent en $O(1)$.

Un million de lectures de 150 pb représente 300 Mo de caractères, sans les objets Python. Une mesure locale avec Python 3.13 donne 191 octets par chaîne de 150 caractères, 272 octets par dictionnaire et 1 256 octets pour la liste des 150 scores, soit 1,9 Ko par lecture matérialisée, près de 2 Go pour un million. Le streaming évite cette facture.

Le fichier gzip réduit la taille sur disque, mais pas la mémoire nécessaire après décompression. La décompression coûte du temps CPU et doit être faite en flux, avec gzip.open(..., "rt"). Les deux octets magiques 1f 8b permettent de reconnaître une archive sans se fier à l'extension.

Implémentation

Les fonctions ci-dessous lisent chaque caractère une fois, vérifient chaque longueur et ne corrigent aucune erreur silencieusement. La version complète et les tests sont dans exercises/module-1/1.2-formats-fasta-fastq/solution_1_2.py ; les sorties sont vérifiées avec Python 3.13.5.

Lire un FASTA

Le parseur conserve les doublons dans une liste, concatène les lignes de séquence et normalise les espaces internes. Il ignore les lignes vides et les commentaires ; avant le premier en-tête. Une séquence avant le premier en-tête, un en-tête vide, un ; à l'intérieur d'un enregistrement ou un > dans une ligne de séquence lèvent une erreur.

def parse_fasta(text):
    records = []
    header = None
    chunks = []
    for line_number, raw_line in enumerate(text.splitlines(), start=1):
        line = raw_line.strip()
        if not line:
            continue
        if line.startswith(">"):
            if header is not None:
                parts = header.split(maxsplit=1)
                records.append({
                    "id": parts[0],
                    "description": parts[1] if len(parts) > 1 else "",
                    "header": header,
                    "sequence": "".join(chunks),
                })
            header = line[1:].strip()
            if not header:
                raise ValueError(f"empty FASTA header at line {line_number}")
            chunks = []
        elif header is None:
            if line.startswith(";"):
                continue
            raise ValueError(f"sequence data before first FASTA header at line {line_number}")
        else:
            if line.startswith(";"):
                raise ValueError(f"legacy ';' comment inside a FASTA record at line {line_number}")
            if ">" in line:
                raise ValueError(f"'>' inside a FASTA sequence line at line {line_number}")
            chunks.append("".join(line.split()))
    if header is not None:
        parts = header.split(maxsplit=1)
        records.append({
            "id": parts[0],
            "description": parts[1] if len(parts) > 1 else "",
            "header": header,
            "sequence": "".join(chunks),
        })
    return records
 
ecoli = """>NC_000913.3:1-60 Escherichia coli str. K-12 substr. MG1655, complete genome
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTC
"""
print(parse_fasta(ecoli))

Sortie vérifiée :

[{'id': 'NC_000913.3:1-60', 'description': 'Escherichia coli str. K-12 substr. MG1655, complete genome', 'header': 'NC_000913.3:1-60 Escherichia coli str. K-12 substr. MG1655, complete genome', 'sequence': 'AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTC'}]

L'identifiant est le premier mot de l'en-tête, la description le reste. La séquence est une liste de morceaux jointe à la fin : sequence += line dans une boucle peut devenir quadratique. La casse et les symboles ambiguus sont conservés : le parseur ne décide pas si la séquence est de l'ADN, de l'ARN ou une protéine.

Lire un FASTQ

Le parseur FASTQ vérifie les quatre lignes, le séparateur +, la longueur de la qualité et l'intervalle ASCII. Il accepte phred33 et phred64. Les lignes vides entre enregistrements sont ignorées ; une ligne vide à l'intérieur d'un enregistrement produit une longueur incohérente ou une lecture vide, jamais une correction silencieuse.

def parse_fastq(text, encoding="phred33"):
    offsets = {"phred33": 33, "phred64": 64}
    if encoding not in offsets:
        raise ValueError(f"unknown encoding: {encoding}")
    offset = offsets[encoding]
    lines = text.splitlines()
    reads = []
    index = 0
    while index < len(lines):
        if not lines[index].strip():
            index += 1
            continue
        if not lines[index].startswith("@"):
            raise ValueError(f"expected '@' header at line {index + 1}")
        header = lines[index][1:].strip()
        if not header:
            raise ValueError(f"empty FASTQ header at line {index + 1}")
        index += 1
        if index + 2 >= len(lines):
            raise ValueError(f"truncated FASTQ record {header.split()[0]!r}")
        sequence = lines[index].strip()
        plus = lines[index + 1]
        quality = lines[index + 2]
        if not plus.startswith("+"):
            raise ValueError(
                f"expected '+' separator for read {header.split()[0]!r} "
                f"at line {index + 2}"
            )
        if len(quality) != len(sequence):
            raise ValueError(
                f"quality length {len(quality)} != sequence length "
                f"{len(sequence)} for read {header.split()[0]!r}"
            )
        scores = [ord(char) - offset for char in quality]
        if any(score < 0 or score > 93 for score in scores):
            raise ValueError(
                f"invalid {encoding} quality character in read "
                f"{header.split()[0]!r}"
            )
        ident, _, description = header.partition(" ")
        reads.append({
            "id": ident,
            "description": description,
            "header": header,
            "sequence": sequence,
            "quality": quality,
            "phred": scores,
        })
        index += 3
    return reads

Sur la lecture read2 du jeu synthétique, on obtient :

{'id': 'read2', 'description': '', 'header': 'read2', 'sequence': 'TTAGGC', 'quality': '##IIII', 'phred': [2, 2, 40, 40, 40, 40]}

La qualité est prise telle quelle, sans strip() : splitlines() a déjà retiré les fins de ligne, et un caractère valide est un ASCII imprimable. La vérification score < 0 or score > 93 couvre les caractères sous l'offset et au-dessus de 126. Pour Phred+64, la borne haute réelle est 62, automatiquement respectée par la limite ASCII.

Le décodage d'une qualité ne demande pas de bibliothèque. La fonction suivante est volontairement minimale ; elle sera réutilisée dans les statistiques et le mini-projet.

Aucune sortie. Exécutez le code pour voir le résultat.

Sortie vérifiée :

Q: [40, 40, 40, 40, 2]
P: [0.0001, 0.0001, 0.0001, 0.0001, 0.630957]
P(no error): 0.368895
P(error): 0.631105
expected errors: 0.631357

Statistiques et streaming

fastq_stats agrège les longueurs et les scores sans supposer des lectures de même longueur. La moyenne de qualité porte sur toutes les bases : une longue lecture de mauvaise qualité pèse plus qu'une courte lecture parfaite. Les fractions Q20/Q30 sont des proportions entre 0 et 1.

def fastq_stats(text, encoding="phred33"):
    reads = parse_fastq(text, encoding=encoding)
    lengths = [len(read["sequence"]) for read in reads]
    total = sum(lengths)
    scores = [score for read in reads for score in read["phred"]]
    q20 = sum(1 for score in scores if score >= 20)
    q30 = sum(1 for score in scores if score >= 30)
    errors = sum(phred_to_probability(score) for score in scores)
    return {
        "n_reads": len(reads),
        "total_bases": total,
        "min_length": min(lengths) if lengths else 0,
        "max_length": max(lengths) if lengths else 0,
        "mean_length": total / len(reads) if reads else 0.0,
        "mean_quality": sum(scores) / total if total else 0.0,
        "q20_bases": q20,
        "q30_bases": q30,
        "q20_fraction": q20 / total if total else 0.0,
        "q30_fraction": q30 / total if total else 0.0,
        "expected_errors": errors,
        "mean_error_probability": errors / total if total else 0.0,
        "encoding": encoding,
    }

Sur le FASTQ synthétique de l'étude de cas, la sortie est :

{'n_reads': 3, 'total_bases': 28, 'min_length': 6, 'max_length': 12, 'mean_length': 9.333333333333334, 'mean_quality': 33.714285714285715, 'q20_bases': 26, 'q30_bases': 26, 'q20_fraction': 0.9285714285714286, 'q30_fraction': 0.9285714285714286, 'expected_errors': 1.2735146889603866, 'mean_error_probability': 0.04548266746287095, 'encoding': 'phred33'}

Pour de vrais fichiers, évitez de charger le texte entier. Le générateur suivant lit un fichier brut ou gzip, détecte le format par les octets magiques et ne stocke que la lecture courante.

from gzip import open as gzip_open
 
def open_text(path):
    with open(path, "rb") as probe:
        magic = probe.read(2)
    if magic == b"\x1f\x8b":
        return gzip_open(path, "rt", encoding="utf-8", newline="")
    return open(path, "rt", encoding="utf-8", newline="")
 
def iter_fastq(path, encoding="phred33"):
    offsets = {"phred33": 33, "phred64": 64}
    offset = offsets[encoding]
    with open_text(path) as handle:
        lines = enumerate(handle, start=1)
        for line_number, raw_line in lines:
            line = raw_line.rstrip("\r\n")
            if not line:
                continue
            if not line.startswith("@"):
                raise ValueError(f"expected '@' at line {line_number}")
            header = line[1:].strip()
            if not header:
                raise ValueError(f"empty FASTA header at line {line_number}")
            try:
                _, sequence_raw = next(lines)
                _, plus_raw = next(lines)
                _, quality_raw = next(lines)
            except StopIteration as exc:
                raise ValueError(f"truncated record {header!r}") from exc
            sequence = sequence_raw.rstrip("\r\n").strip()
            plus = plus_raw.rstrip("\r\n")
            quality = quality_raw.rstrip("\r\n")
            if not plus.startswith("+"):
                raise ValueError(f"missing '+' in record {header!r}")
            if len(quality) != len(sequence):
                raise ValueError(
                    f"quality length {len(quality)} != sequence length "
                    f"{len(sequence)} for {header!r}"
                )
            scores = [ord(char) - offset for char in quality]
            if any(score < 0 or score > 93 for score in scores):
                raise ValueError(f"invalid quality character for {encoding}")
            ident, _, description = header.partition(" ")
            yield {
                "id": ident,
                "description": description,
                "header": header,
                "sequence": sequence,
                "quality": quality,
                "phred": scores,
            }

Le générateur donne les mêmes identifiants pour le fichier brut et pour sa version gzip :

plain: ['read1', 'read2', 'read3']
gzip: ['read1', 'read2', 'read3']

newline="" empêche la traduction des fins de ligne et laisse rstrip("\r\n") faire le travail. splitlines() est robuste, mais charge tout le texte ; le générateur garde la mémoire bornée par la plus grande lecture.

Écrire et vérifier l'aller-retour

Écrire un FASTQ est simple si l'on conserve la qualité exacte. Pour un FASTA, le repli de ligne est une présentation : un aller-retour compare les enregistrements, pas les octets.

def format_fasta(records, width=60):
    lines = []
    for record in records:
        lines.append(">" + record["header"])
        sequence = record.get("sequence", "")
        for start in range(0, len(sequence), width):
            lines.append(sequence[start:start + width])
    return "\n".join(lines) + ("\n" if lines else "")
 
 
def format_fastq(records):
    lines = []
    for record in records:
        sequence = record["sequence"]
        quality = record["quality"]
        if len(sequence) != len(quality):
            raise ValueError(
                f"quality length {len(quality)} != sequence length "
                f"{len(sequence)} for read {record.get('id', '?')!r}"
            )
        lines.extend(["@" + record["header"], sequence, "+", quality])
    return "\n".join(lines) + ("\n" if lines else "")
 
 
fastq_records = [
    {"header": "read1", "sequence": "ACGT", "quality": "IIII"},
    {"header": "read2", "sequence": "TT", "quality": "??"},
]
fastq_text = format_fastq(fastq_records)
print(repr(fastq_text))
print(fastq_text == "@read1\nACGT\n+\nIIII\n@read2\nTT\n+\n??\n")
 
fasta_records = [{"header": "seq1", "sequence": "ACGTAC"}]
print(repr(format_fasta(fasta_records, width=4)))

Sortie vérifiée :

'@read1\nACGT\n+\nIIII\n@read2\nTT\n+\n??\n'
True
'>seq1\nACGT\nAC\n'

Le même test appliqué aux enregistrements issus de parse_fastq donne True, et parse_fasta(format_fasta(records)) == records est également True : le repli de ligne et la ligne + sont normalisés, mais les champs utiles sont conservés.

Écrivez dans un fichier temporaire puis renommez, pour éviter un fichier tronqué. Pour gzip, utilisez newline="\n" afin d'obtenir des fins de ligne reproductibles sur Windows comme sur Linux.

Les erreurs de fichiers réels

Les défauts fréquents sont la longueur de qualité incohérente, le séparateur + absent, l'enregistrement tronqué et les fins de ligne CRLF. split("\n") garde un \r et ajoute un caractère à la qualité ; splitlines() retire correctement \r\n.

crlf = "@r1\r\nACGT\r\n+\r\nIIII\r\n"
print("split:", [len(line) for line in crlf.split("\n")])
print("splitlines:", [len(line) for line in crlf.splitlines()])
print("ids:", [read["id"] for read in parse_fastq(crlf)])
 
for label, text in {
    "short quality": "@r1\nACGT\n+\nIII\n",
    "missing plus": "@r1\nACGT\n-\nIIII\n",
    "truncated": "@r1\nACGT\n+\n",
}.items():
    try:
        parse_fastq(text)
    except ValueError as exc:
        print(label, "->", exc)

Sortie vérifiée :

split: [4, 5, 2, 5, 0]
splitlines: [3, 4, 1, 4]
ids: ['r1']
short quality -> quality length 3 != sequence length 4 for read 'r1'
missing plus -> expected '+' separator for read 'r1' at line 3
truncated -> truncated FASTQ record 'r1'

Un en-tête dupliqué n'est pas une erreur de syntaxe. Convertir la liste en dictionnaire indexé par id écrase le premier ; le parseur renvoie donc une liste, et le mini-projet compte les doublons.

Étude de cas

Nous travaillons sur deux jeux : un FASTQ synthétique de trois lectures, étiqueté comme exemple, et un extrait FASTA réel du génome d'E. coli K-12 MG1655. Le premier permet de contrôler les longueurs et la qualité ; le second montre ce qu'un FASTA peut et ne peut pas dire.

@read1
ACGTACGTACGT
+
IIIIIIIIIIII
@read2
TTAGGC
+
##IIII
@read3
GGGCCCAAAT
+
??????????

La question est double : les longueurs de qualité sont-elles cohérentes, et les lectures passent-elles un contrôle raisonnable ? Nous appliquons parse_fastq, puis fastq_stats.

Les trois longueurs de qualité valent 12, 6 et 10, exactement les longueurs des séquences : aucun défaut structurel. Les statistiques donnent 28 bases, 26 bases à Q20 ou plus, soit 92,9 %, et une qualité moyenne de 33,7. Ces chiffres agrégés sont bons. Pourtant, read2 contient deux bases à Q2 : sa probabilité d'au moins une erreur vaut $1-(1-0{,}630957)^2 \approx 0{,}864$, soit 86 %. La moyenne de 33,7 masque ce défaut, car les 12 bases parfaites de read1 et les 10 bases Q30 de read3 compensent les deux bases catastrophiques.

La décision dépend du seuil fixé à l'avance. Si le protocole exige toutes les bases à Q20, read2 est rejetée ou trimée ; sinon, elle peut être conservée avec un avertissement. Le nombre attendu d'erreurs, 1,27 sur 28 bases, confirme que le jeu n'est pas aussi propre que la moyenne ne le suggère.

L'extrait FASTA a été obtenu avec la commande suivante, exécutée le 2026-11-01 :

efetch -db nuccore -id NC_000913.3 -seq_start 1 -seq_stop 60 -format fasta

Le parseur renvoie un identifiant NC_000913.3:1-60, une description complète et une séquence de longueur 60. Les comptes sont A = 17, C = 10, G = 14, T = 19, soit un GC de 40,0 %. Ce chiffre est inférieur au GC global d'E. coli K-12, de l'ordre de 50 %, mais une fenêtre de 60 pb est trop courte pour conclure : cette région peut être localement AT-riche.

Le FASTA ne contient aucune qualité : on peut vérifier la longueur, l'alphabet et la composition, pas la probabilité d'erreur de chaque base. Cette limite est structurelle, pas un défaut du parseur. Pour contrôler la qualité, il faudrait le FASTQ dont cet extrait est issu, ou un fichier QUAL associé. La leçon méthodologique est de ne pas demander à un format une information qu'il ne stocke pas.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-1/1.2-formats-fasta-fastq -q

python -m pytest exercises/module-1/1.2-formats-fasta-fastq -q

python -m pytest exercises/module-1/1.2-formats-fasta-fastq/projet -q

Résumé

• Un FASTA contient un en-tête >, un identifiant (premier mot) et une ou plusieurs lignes de séquence ; il ne contient pas de qualité.
• Un FASTQ moderne tient sur quatre lignes : @, séquence, +, qualité de même longueur que la séquence.
• Le score Phred vérifie $Q = -10\log_{10}P$ ; Q20 correspond à 1 % d'erreur, Q30 à 0,1 %.
• La moyenne des $Q$ n'est pas une probabilité : une seule base à Q2 peut rendre une lecture peu fiable.
• Phred+33 est l'encodage moderne ; Phred+64 est historique et change la valeur de I.
• Un parseur robuste conserve les doublons, gère CRLF et les lignes vides, vérifie les longueurs et lit gzip en streaming.
• Le coût est $O(L)$ en temps et $O(L)$ en mémoire si l'on matérialise tout ; le streaming réduit la mémoire à la taille d'une lecture.
• Un aller-retour lecture-écriture se vérifie en comparant les enregistrements, pas les octets.

Références

• Cock, P. J. A. et al. (2009). Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics, 25(11), 1422–1423.
• Chen, S., Zhou, Y., Chen, Y. & Gu, J. (2018). fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics, 34(17), i884–i890.
• NCBI. FASTA Format for Nucleotide Sequences. https://www.ncbi.nlm.nih.gov/genbank/fastaformat/
• Biopython. Bio.SeqIO package (documentation 1.88). https://biopython.org/docs/latest/api/Bio.SeqIO.html
• NCBI RefSeq. Escherichia coli str. K-12 substr. MG1655, complete genome, accession NC_000913.3. https://www.ncbi.nlm.nih.gov/nuccore/NC_000913.3
• Python. gzip — Support for gzip files. https://docs.python.org/3/library/gzip.html