Parsing avec Biopython SeqIO

Prérequis

Cette section prolonge 1.1 et surtout 1.2 : vous devez reconnaître un en-tête FASTA, une ligne de séquence et les quatre lignes d'un FASTQ. Côté Python, fonctions, dictionnaires, compréhensions et pathlib doivent être familiers. Biopython 1.88 est supposé installé ; le code est compatible avec Python 3.11 et versions ultérieures.

Introduction

Un séquenceur produit des fichiers, pas une séquence abstraite. Le génome de référence d'Escherichia coli K-12 MG1655, accession NC_000913.3, compte environ 4,6 millions de paires de bases ; une expérience d'ARN-seq peut produire des dizaines de millions de lectures, chacune avec un score de qualité par base. Entre le fichier brut et l'analyse, il faut lire, valider, filtrer, convertir et réécrire. Une séquence tronquée, un identifiant dupliqué ou une qualité mal interprétée fausse silencieusement les comptages en aval.

Le défi est plus subtil qu'il n'y paraît. Un FASTA autorise plusieurs lignes par séquence, des descriptions après l'identifiant, des lignes vides et des fins de ligne différentes. Un FASTQ ajoute une ligne de qualité dont la longueur doit correspondre exactement à la séquence, avec un encodage ASCII qui a changé selon les plateformes. Un parseur maison devient risqué dès qu'il faut gérer ces variantes, les formats annotés et les fichiers volumineux.

Biopython fournit Bio.SeqIO, une interface unique pour une trentaine de formats. Cette section montre quand appeler parse plutôt que read, ce que contient un SeqRecord, pourquoi record.seq n'est pas une chaîne ordinaire, comment écrire FASTA et FASTQ, et comment indexer un fichier trop gros pour être chargé. Le fil rouge reste NC_000913.3, avec un extrait réel de 60 pb et un FASTQ synthétique.

Théorie

Du texte brut à l'objet structuré

Un FASTA semble se réduire à deux types de lignes : un en-tête commençant par > et des lignes de séquence. L'identifiant est le premier mot de l'en-tête, mais la description peut contenir espaces, crochets et numéros de version. Une séquence peut être coupée en lignes de 60, 70 ou 80 caractères, et des minuscules peuvent signaler un masquage. Un split(">") naïf produit donc des objets faux dès que le fichier s'écarte du cas d'école.

Un parseur implémente une grammaire : il reconnaît les en-têtes, concatène les lignes, associe les qualités aux bases et conserve les métadonnées. SeqIO fait ce travail pour FASTA, FASTQ, GenBank, EMBL, UniProt et d'autres formats. La même boucle fonctionne sur un FASTA et sur un GenBank.

parse, read et les autres modes d'accès

SeqIO.parse(handle, format) renvoie un itérateur. Il ne lit pas tout d'un coup : chaque next produit le SeqRecord suivant. La mémoire reste proportionnelle au plus gros enregistrement, pas au fichier. En contrepartie, l'itérateur ne se parcourt qu'une fois et les erreurs de format peuvent survenir pendant l'itération.

SeqIO.read(handle, format) renvoie exactement un SeqRecord. Avec zéro ou plusieurs enregistrements, il lève ValueError. Il convient à un plasmide ou à un petit GenBank. read ne signifie donc pas « lire tout », mais « lire un seul enregistrement, et exiger qu'il soit seul ».

SeqIO.to_dict charge tous les enregistrements dans un dictionnaire indexé par identifiant. Pratique sur un petit fichier, coûteux en mémoire, et interdit les doublons : deux en-têtes identiques produisent ValueError: Duplicate key. Le tableau 1.3.1 résume les modes d'accès.

La complexité est linéaire en caractères pour parse, read et to_dict. La construction d'un index est aussi linéaire, mais une consultation par identifiant devient quasi constante.

Anatomie d'un SeqRecord

Un SeqRecord est le conteneur central. Ses attributs principaux sont id, name, description, seq, features, annotations, letter_annotations et dbxrefs. Pour un FASTA, id vaut le premier mot de l'en-tête, name reprend souvent cet identifiant, et description contient l'en-tête complet. Pour un GenBank, les annotations et les features décrivent gènes, CDS et promoteurs.

Un SeqFeature possède un type, une location et des qualifiers. Dans Biopython 1.88, le brin se lit sur feature.location.strand, et non sur feature.strand. Les coordonnées location.start et location.end sont 0-based et semi-ouvertes : [start, end) couvre start à end - 1. Pour un rapport 1-based inclusif, affichez start + 1 et end. feature.extract(record.seq) renvoie la sous-séquence et applique le complément inverse si le brin vaut -1.

Un FASTA ne contient aucune feature : record.features est vide. Pour des annotations, téléchargez un GenBank ou croisez les coordonnées avec un GTF/GFF ; aucun parseur ne peut inventer une annotation absente.

Seq n'est pas str

record.seq est un Bio.Seq.Seq, pas une chaîne. Il partage len, l'indexation, le découpage, count et find, mais un découpage renvoie un Seq, et l'objet possède reverse_complement, complement, transcribe et translate. Dans 1.88, Seq("ACGT") == "ACGT" vaut True : l'égalité ne suffit pas à les distinguer.

La différence pratique est ailleurs : isinstance(record.seq, str) est faux, json.dumps refuse un Seq, et certaines bibliothèques compilées rejettent l'objet. Convertissez avec str(record.seq) aux frontières : JSON, fichier texte, expression régulière, hash. À l'intérieur, gardez le Seq pour ses méthodes et éviter les copies. Seq est immuable ; MutableSeq sert aux modifications en place.

Écrire, convertir et préserver les qualités

SeqIO.write(records, handle, format) écrit un itérable et renvoie le nombre d'objets écrits. Le format vaut "fasta", "fastq" ou "genbank". La fonction accepte un générateur, ce qui évite de matérialiser la liste complète ; elle ne ferme pas un handle ouvert par l'appelant.

Pour contrôler le FASTA, FastaIO.FastaWriter(handle, wrap=70) expose write_record et write_file. wrap règle la largeur des lignes ; la valeur par défaut est 60. Pour le FASTQ, Bio.SeqIO.FastqIO n'existe pas en 1.88 : l'import lève ModuleNotFoundError. L'écriture de bas niveau est dans Bio.SeqIO.QualityIO.FastqPhredWriter, la lecture détaillée dans FastqGeneralIterator.

Un enregistrement FASTQ exige un score par base dans record.letter_annotations["phred_quality"]. Sans cette liste, SeqIO.write(..., "fastq") lève ValueError: No suitable quality scores found.... La conversion FASTA vers FASTQ est donc impossible : le FASTA ne porte aucune qualité. Inventer des qualités produit un fichier syntaxiquement valide mais scientifiquement faux. La conversion inverse, FASTQ vers FASTA, est légitime mais perd les qualités, ce qui doit être documenté.

Le score Phred $Q$ relie la probabilité d'erreur $P$ à une échelle logarithmique : $Q = -10\log_{10} P$. En encodage Sanger et Illumina 1.8+, le caractère ASCII vaut $Q + 33$ ; I code 40, ? code 30 et + code 10. Un score de 30 correspond à $10^{-3}$.

GC, masse moléculaire et cadres de lecture

Bio.SeqUtils.gc_fraction(seq) renvoie une fraction entre 0 et 1, et non un pourcentage. Pour une séquence $s$, la définition est

Le paramètre ambiguous change le dénominateur. La valeur par défaut "remove" retire les symboles ambigus ; "ignore" les garde dans le dénominateur sans les compter comme G ou C ; "weighted" leur donne une contribution fractionnaire. gc_fraction("ACGTN") vaut 0,5 par défaut et 0,4 avec "ignore". Une séquence entièrement ambiguë renvoie 0,0. Précisez toujours la convention dans un rapport.

molecular_weight(seq, seq_type="DNA") calcule une masse en daltons. Par défaut, l'ADN est simple brin ; double_stranded=True donne le duplex, circular=True retire la molécule d'eau de fermeture, et seq_type="protein" traite une protéine. La masse est une somme de masses de résidus corrigée des condensations ; elle ne remplace pas une mesure expérimentale. La fonction lève ValueError sur un symbole non ambigu comme N.

six_frame_translations(seq) produit un rapport textuel des six cadres, avec comptage des bases et GC. C'est un outil d'inspection, pas une structure de données : dans un pipeline, traduisez les six séquences avec Seq.translate et conservez les coordonnées. Le rapport reste utile pour repérer un cadre ouvert ou un codon stop.

Exemple travaillé à la main

Reprenons l'extrait réel de 60 pb au début de NC_000913.3 :

AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTC

Le comptage donne A = 17, C = 10, G = 14 et T = 19. Le GC vaut $(14 + 10) / 60 = 0{,}4$, soit 40 %. La masse du simple brin est 18 622,9 Da, celle du duplex 37 103,72 Da, soit environ 310 Da par nucléotide et 618 Da par paire de bases.

Le second enregistrement, synthétique, fait 39 pb : ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG. A = 9, C = 8, G = 14, T = 8, donc un GC de $22/39 \approx 56{,}41\,\%$ et une masse de 12 192,76 Da. Ces valeurs serviront de contrôle.

Indexer plutôt que tout charger

SeqIO.index(filename, format) indexe les décalages du fichier et renvoie un objet semblable à un dictionnaire. Il ne charge pas les séquences : index["id"] parse seulement l'enregistrement demandé. Il exige des identifiants uniques et un nom de fichier, pas un handle. index_db(db_path, filenames, format) stocke l'index dans SQLite et le réutilise entre deux exécutions. Les deux conviennent aux fichiers de plusieurs gigaoctets ; to_dict non.

Un index devient faux si le fichier change : réindexez après modification et fermez l'objet avec index.close(). Pour un parcours complet, parse reste plus simple ; l'index n'est utile que si les identifiants demandés sont connus.

Implémentation

Vérifions d'abord la version installée. Les sorties de cette section ont été exécutées le 2026-11-01 avec Biopython 1.88 et Python 3.13.

python -c "import Bio; print(Bio.__version__)"

1.88

Lire deux enregistrements et calculer GC et masse

Le bloc suivant construit un FASTA en mémoire, puis parcourt ses deux enregistrements. Le premier est l'extrait réel de NC_000913.3 ; le second est un exemple synthétique, clairement étiqueté.

from io import StringIO
 
from Bio import SeqIO
from Bio.SeqUtils import gc_fraction, molecular_weight
 
fasta_text = """>NC_000913.3:1-60 E. coli K-12 MG1655, positions 1-60
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTC
>synth_39 synthetic example, not a real locus
ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG
"""
 
for record in SeqIO.parse(StringIO(fasta_text), "fasta"):
    gc = 100 * gc_fraction(record.seq)
    mass = molecular_weight(record.seq, seq_type="DNA")
    print(f"{record.id:22s} {len(record.seq):3d} pb  GC {gc:5.1f} %  {mass:10.1f} Da")

NC_000913.3:1-60        60 pb  GC  40.0 %     18622.9 Da
synth_39                39 pb  GC  56.4 %     12192.8 Da

Les valeurs confirment le calcul manuel. gc_fraction accepte directement un Seq ; le formatage :22s aligne les identifiants.

Itérateur, lecture unique et erreurs

SeqIO.parse renvoie un itérateur. La preuve la plus simple est de l'épuiser puis de recommencer la boucle.

from io import StringIO
 
from Bio import SeqIO
 
handle = StringIO(">a\nACGT\n>b\nTTTT\n")
records = SeqIO.parse(handle, "fasta")
print(type(records).__name__)
print(next(records).id)
print([record.id for record in records])
print([record.id for record in records])
 
try:
    SeqIO.read(StringIO(">a\nACGT\n>b\nTTTT\n"), "fasta")
except ValueError as error:
    print(type(error).__name__, ":", error)

FastaIterator
a
['b']
[]
ValueError : More than one record found in handle

Le deuxième parcours ne renvoie rien : l'itérateur est à usage unique. Pour deux parcours, transformez-le en liste avec list(...). SeqIO.read échoue sur deux enregistrements ; pour un fichier inconnu, préférez parse et comptez les objets.

Seq, str et sérialisation

Le bloc suivant montre les opérations qui fonctionnent sur un Seq et le point où la conversion devient nécessaire.

from json import dumps
 
from Bio.Seq import Seq
 
seq = Seq("AGCTTTTCATT")
print(type(seq).__name__, len(seq))
print(type(seq[:4]).__name__, seq[:4])
print(seq.reverse_complement())
print(str(seq))
print(seq == "AGCTTTTCATT")
 
try:
    dumps({"sequence": seq})
except TypeError as error:
    print(type(error).__name__, ":", error)

Seq 11
Seq AGCT
AATGAAAAGCT
AGCTTTTCATT
True
TypeError : Object of type Seq is not JSON serializable

Le découpage renvoie encore un Seq et l'égalité avec une chaîne vaut True, mais JSON échoue. Gardez le Seq pour les méthodes biologiques et convertissez avec str() dès que l'objet sort du programme.

Écrire un FASTA, puis buter sur le FASTQ

Le bloc ci-dessous écrit deux enregistrements en FASTA et tente la conversion en FASTQ. La seconde écriture échoue, et c'est précisément le comportement à comprendre.

from pathlib import Path
from tempfile import TemporaryDirectory
 
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
 
records = [
    SeqRecord(
        Seq("AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTC"),
        id="NC_000913.3:1-60",
        description="E. coli K-12 MG1655, positions 1-60",
    ),
    SeqRecord(
        Seq("ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG"),
        id="synth_39",
        description="synthetic example",
    ),
]
 
with TemporaryDirectory() as tmp:
    fasta_path = Path(tmp) / "records.fasta"
    print("written:", SeqIO.write(records, fasta_path, "fasta"))
    print(fasta_path.read_text(encoding="utf-8"))
    try:
        SeqIO.write(records, Path(tmp) / "records.fastq", "fastq")
    except ValueError as error:
        print(type(error).__name__, ":", error)

written: 2
>NC_000913.3:1-60 E. coli K-12 MG1655, positions 1-60
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTC
>synth_39 synthetic example
ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG
 
ValueError : No suitable quality scores found in letter_annotations of SeqRecord (id=NC_000913.3:1-60).

Le FASTA est écrit avec les descriptions complètes. La tentative FASTQ échoue avant d'écrire : Biopython refuse de fabriquer des qualités absentes. Pour un FASTQ, ajoutez letter_annotations["phred_quality"].

from pathlib import Path
from tempfile import TemporaryDirectory
 
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
 
record = SeqRecord(
    Seq("AGCTTTTCATTCTGACTGCAACGGGCAA"),
    id="read_001",
    description="synthetic read",
)
record.letter_annotations["phred_quality"] = [40] * len(record.seq)
 
with TemporaryDirectory() as tmp:
    path = Path(tmp) / "read.fastq"
    SeqIO.write([record], path, "fastq")
    print(path.read_text(encoding="utf-8"))
    back = SeqIO.read(path, "fastq")
    print(back.id, back.letter_annotations["phred_quality"][:5], len(back))

@read_001 synthetic read
AGCTTTTCATTCTGACTGCAACGGGCAA
+
IIIIIIIIIIIIIIIIIIIIIIIIIIII
 
read_001 [40, 40, 40, 40, 40] 28

I code bien 40 et la relecture restitue des entiers, pas des caractères ASCII. Pour un contrôle fin du FASTA, FastaIO.FastaWriter accepte wrap=70. Pour le FASTQ de bas niveau, utilisez QualityIO.FastqPhredWriter ; Bio.SeqIO.FastqIO n'existe pas.

try:
    from Bio.SeqIO import FastqIO
except ModuleNotFoundError as error:
    print(type(error).__name__, ":", error)
 
from Bio.SeqIO.QualityIO import FastqPhredWriter
 
print(FastqPhredWriter.__name__)

ModuleNotFoundError : No module named 'Bio.SeqIO.FastqIO'
FastqPhredWriter

GC, masse et six cadres

Le bloc suivant rassemble les paramètres qui changent les résultats numériques.

from Bio.SeqUtils import gc_fraction, molecular_weight, six_frame_translations
 
seq = "AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTC"
 
print(gc_fraction(seq))
print(round(100 * gc_fraction(seq), 1))
print(round(molecular_weight(seq), 1))
print(round(molecular_weight(seq, double_stranded=True), 1))
print(gc_fraction("ACGTN"))
print(gc_fraction("ACGTN", ambiguous="ignore"))
print(gc_fraction("NNNN"))

0.4
40.0
18622.9
37103.7
0.5
0.4
0.0

La différence entre 0,5 et 0,4 pour ACGTN n'est pas un bug : elle vient du dénominateur. Le mode par défaut retire le N, le mode "ignore" le compte dans la longueur totale. Dans un rapport, écrivez toujours « GC calculé après retrait des bases ambiguës » ou « GC sur la longueur totale ».

Le rapport de six cadres est un texte, utile pour l'inspection visuelle. En voici la sortie réelle pour l'extrait de 60 pb.

from Bio.SeqUtils import six_frame_translations
 
seq = "AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTC"
print(six_frame_translations(seq))

GC_Frame: a:17 t:19 g:14 c:10
Sequence: agcttttcat ... aaagagtgtc, 60 nt, 40.00 %GC
 
 
1/1
  L  F  I  L  T  A  T  G  N  M  S  L  C  G  L  K  K  E  C
 A  F  H  S  D  C  N  G  Q  Y  V  S  V  W  I  K  K  R  V
S  F  S  F  *  L  Q  R  A  I  C  L  C  V  D  *  K  K  S  V
agcttttcattctgactgcaacgggcaatatgtctctgtgtggattaaaaaaagagtgtc   40 %
tcgaaaagtaagactgacgttgcccgttatacagagacacacctaatttttttctcacag
K  E  N  Q  S  C  R  A  I  H  R  Q  T  S  *  F  F  L  T
 A  K  *  E  S  Q  L  P  C  Y  T  E  T  H  I  L  F  L  T  D
  S  K  M  R  V  A  V  P  L  I  D  R  H  P  N  F  F  S  H

Le rapport affiche le comptage, la séquence en minuscules, puis les cadres directs et inverses. Ce n'est pas une prédiction de gène : un cadre ouvert peut être fortuit.

Indexer un fichier

Pour un fichier que l'on interroge par identifiant, SeqIO.index évite de charger toutes les séquences. L'exemple crée un petit FASTA, l'indexe, lit un enregistrement par son identifiant, puis construit un index SQLite persistant.

from pathlib import Path
from tempfile import TemporaryDirectory
 
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
 
records = [
    SeqRecord(Seq("AGCTTTTCATTCTGACTGCAACGGGCAA"), id="ecoli_1_60"),
    SeqRecord(Seq("ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG"), id="synth_39"),
]
 
with TemporaryDirectory(ignore_cleanup_errors=True) as tmp:
    fasta_path = Path(tmp) / "records.fasta"
    SeqIO.write(records, fasta_path, "fasta")
 
    index = SeqIO.index(fasta_path, "fasta")
    print(len(index), list(index))
    print(index["synth_39"].id, len(index["synth_39"].seq))
    print(index.get_raw("synth_39").decode().splitlines()[0])
    index.close()
 
    db_path = Path(tmp) / "records.sqlite"
    SeqIO.index_db(db_path, fasta_path, "fasta").close()
    reopened = SeqIO.index_db(db_path, fasta_path, "fasta")
    print(list(reopened))
    reopened.close()

2 ['ecoli_1_60', 'synth_39']
synth_39 39
>synth_39 <unknown description>
['ecoli_1_60', 'synth_39']

get_raw renvoie les octets bruts, sans parser la séquence. Notez <unknown description> : un SeqRecord sans description reçoit cette valeur, qui se retrouve dans le FASTA ; passez description="" pour un en-tête propre.

Diagnostiquer les erreurs

Un fichier vide ou multi-enregistrements passé à read produit les messages suivants.

from io import StringIO
 
from Bio import SeqIO
 
for text in ["", ">a\nACGT\n>b\nTTTT\n"]:
    try:
        SeqIO.read(StringIO(text), "fasta")
    except ValueError as error:
        print(type(error).__name__, ":", error)
 
try:
    list(SeqIO.parse(StringIO(">a\nACGT\n"), "xyz"))
except ValueError as error:
    print(type(error).__name__, ":", error)

ValueError : No records found in handle
ValueError : More than one record found in handle
ValueError : Unknown format 'xyz'

Un FASTQ tronqué produit ValueError: Lengths of sequence and quality values differs..., un fichier absent FileNotFoundError, et molecular_weight sur N une ValueError. Ces messages sont des contrats : testez-les plutôt que de les masquer par un except Exception silencieux.

Étude de cas

Nous travaillons sur deux jeux complémentaires. Le premier est un extrait réel du chromosome d'E. coli K-12 MG1655, accession NC_000913.3, positions 1 à 60. Le second est un FASTQ synthétique construit à partir de cet extrait et d'une lecture courte ; les séquences sont réelles ou explicitement synthétiques, mais les scores de qualité du FASTQ sont inventés. Aucune conclusion biologique ne doit être tirée de ces qualités.

La question est la suivante : peut-on lire cet extrait, calculer ses propriétés, le convertir en FASTQ, puis vérifier que la relecture restitue exactement l'information ? La réponse attendue est oui pour le FASTA et non pour une conversion FASTA vers FASTQ sans qualités.

La région a été téléchargée avec EFetch :

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

Le fichier obtenu contient l'en-tête et la séquence suivants.

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

Première étape : lire le FASTA et vérifier les propriétés calculées. Le script ci-dessous est exécuté sur le fichier téléchargé ; les valeurs sont réelles.

from Bio import SeqIO
from Bio.SeqUtils import gc_fraction, molecular_weight
 
record = SeqIO.read("ecoli_1_60.fasta", "fasta")
print(record.id)
print(record.description)
print(len(record.seq))
print(round(100 * gc_fraction(record.seq), 2))
print(round(molecular_weight(record.seq, double_stranded=True), 2))

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

Le GC de 40 % est inférieur à la moyenne du génome, un peu au-dessus de 50 %, mais un window de 60 pb est trop court pour conclure quoi que ce soit sur une région. La masse du duplex, 37 103,72 Da, correspond à l'ordre de grandeur attendu pour 60 pb.

Deuxième étape : construire un FASTQ synthétique. La première lecture reprend l'extrait de 60 pb avec une qualité constante de 40. La seconde est une lecture courte de 30 pb avec des qualités variables, de 40 à 10, pour rendre la lecture des scores visible.

@ecoli_1_60
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTC
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
@synth_lowq
ATGAAACGCATTAGCACCACCATTACCACC
+
IIIIIIIIII??????5555555555++++

Le fichier est ensuite parsé comme n'importe quel FASTQ. Les qualités sont disponibles dans letter_annotations.

from io import StringIO
 
from Bio import SeqIO
 
fastq_text = """@ecoli_1_60
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTC
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
@synth_lowq
ATGAAACGCATTAGCACCACCATTACCACC
+
IIIIIIIIII??????5555555555++++
"""
 
for record in SeqIO.parse(StringIO(fastq_text), "fastq"):
    qualities = record.letter_annotations["phred_quality"]
    mean = sum(qualities) / len(qualities)
    print(record.id, len(record.seq), min(qualities), round(mean, 1))

ecoli_1_60 60 40 40.0
synth_lowq 30 10 27.3

Troisième étape : tenter la conversion FASTA vers FASTQ. Le message d'erreur est explicite et intervient avant toute écriture.

from Bio import SeqIO
 
records = list(SeqIO.parse("ecoli_1_60.fasta", "fasta"))
try:
    SeqIO.write(records, "ecoli_1_60.fastq", "fastq")
except ValueError as error:
    print(type(error).__name__, ":", error)

ValueError : No suitable quality scores found in letter_annotations of SeqRecord (id=NC_000913.3:1-60).

Quatrième étape : écrire le FASTQ à partir des enregistrements qui possèdent des qualités, puis relire et comparer. La comparaison porte sur l'identifiant, la séquence et la liste de qualités.

from io import StringIO
 
from Bio import SeqIO
 
fastq_text = """@ecoli_1_60
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTC
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
@synth_lowq
ATGAAACGCATTAGCACCACCATTACCACC
+
IIIIIIIIII??????5555555555++++
"""
 
records = list(SeqIO.parse(StringIO(fastq_text), "fastq"))
buffer = StringIO()
written = SeqIO.write(records, buffer, "fastq")
buffer.seek(0)
back = list(SeqIO.parse(buffer, "fastq"))
same = all(
    (before.id, str(before.seq), before.letter_annotations["phred_quality"])
    == (after.id, str(after.seq), after.letter_annotations["phred_quality"])
    for before, after in zip(records, back)
)
print(written, len(back), same)

2 2 True

La relecture est identique. Ce test de round-trip ne prouve pas que les qualités sont biologiquement valides ; il prouve que la sérialisation et la désérialisation conservent l'information. C'est exactement le niveau de garantie qu'un pipeline doit vérifier avant de supprimer le fichier intermédiaire.

Cinquième étape : indexer le FASTA et retrouver un enregistrement par son identifiant. Sur un fichier de 60 pb, l'index est superflu, mais la syntaxe est la même que pour un génome complet.

from Bio import SeqIO
 
index = SeqIO.index("ecoli_1_60.fasta", "fasta")
print(list(index))
print(index["NC_000913.3:1-60"].id)
print(len(index["NC_000913.3:1-60"].seq))
index.close()

['NC_000913.3:1-60']
NC_000913.3:1-60
60

L'interprétation biologique reste limitée : l'extrait est réel, mais 60 pb ne suffisent pas à caractériser une région. L'intérêt de l'étude est méthodologique. Un FASTA et un FASTQ ne portent pas la même information ; une conversion qui « marche » syntaxiquement peut avoir perdu ou inventé des qualités. Le contrôle de round-trip, la vérification des identifiants et la documentation des conventions de GC sont les trois garde-fous à conserver dans un pipeline.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-1/1.3-biopython-seqio -q

python -m pytest exercises/module-1/1.3-biopython-seqio -q

python -m pytest exercises/module-1/1.3-biopython-seqio/projet -q

Résumé

• SeqIO.parse renvoie un itérateur paresseux, à usage unique ; SeqIO.read exige exactement un enregistrement.
• Un SeqRecord conserve l'identifiant, la description, la séquence, les features, les annotations et les qualités par base.
• Seq n'est pas une chaîne : str(record.seq) est nécessaire pour JSON, les fichiers texte et les bibliothèques qui exigent strictement str.
• SeqIO.write écrit un itérable et renvoie le nombre d'enregistrements ; FastaIO règle le retour à la ligne, QualityIO gère le FASTQ.
• Un FASTA ne peut pas devenir un FASTQ sans scores de qualité ; la conversion inverse est possible mais perd les qualités.
• gc_fraction renvoie une fraction et traite les bases ambiguës selon le paramètre ambiguous ; molecular_weight refuse les symboles non ambigus.
• SeqIO.index et index_db évitent de charger tout le fichier, mais exigent des identifiants uniques et une réindexation après modification.
• Un test de round-trip vérifie la fidélité de la sérialisation, pas la validité biologique des données.

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.
• Cock, P. J. A., Fields, C. J., Goto, N., Heuer, M. L. et Rice, P. M. (2010). The Sanger FASTQ file format for sequences with quality scores, and the Solexa/Illumina FASTQ variants. Nucleic Acids Research, 38(6), 1767–1771.
• Biopython. Tutorial and Cookbook et documentation de Bio.SeqIO. https://biopython.org/docs/latest/Tutorial/ et https://biopython.org/docs/latest/api/Bio.SeqIO.html
• Biopython. Bio.SeqUtils : gc_fraction, molecular_weight, six_frame_translations. https://biopython.org/docs/latest/api/Bio.SeqUtils.html
• NCBI. Escherichia coli str. K-12 substr. MG1655, complete genome, accession NC_000913.3. https://www.ncbi.nlm.nih.gov/nuccore/NC_000913.3
• NCBI. E-utilities et format FASTA. https://www.ncbi.nlm.nih.gov/books/NBK25501/ et https://www.ncbi.nlm.nih.gov/BLAST/fasta.shtml