> Modules non standards > Biopython > Lecture et écriture des séquences
Lecture et écriture des séquences
Lecture de séquences :
- import : from Bio import SeqIO
- lecture d'un fichier fasta : for seqRecord in SeqIO.parse('myFile.fas', 'fasta') ('genbank' pour 'genbank') : renvoie des objets SeqRecord
- on peut faire la même chose avec un fichier fastq : for seqRecord in SeqIO.parse('myFile.fastq', 'fastq'), mais attention, ce n'est pas ce qu'il y a de plus rapide (5 fois plus lent que lecture manuelle).
- on peut aussi donner un filehandle : SeqIO.parse(myFh, 'genbank')
- beaucoup de formats supportés : fasta, genbank, embl, clustal, fastaq, sff, swiss (uniprot), qual (qualité en format fasta), uniprot-xml, tab (simple tableau) ...
- si le fichier contient une seule séquence, on peut faire simplement : seqRecord = SeqIO.read('myFile.fas', 'fasta')
- on peut faire aussi : it = SeqIO.parse('myFile.gbk', 'genbank'); seq1 = next(it), seq2 = next(it); .... Il faut alors faire :
it = SeqIO.parse('toto.fastq', 'fastq')
while True:
try:
seqRecord = next(it)
print(seqRecord.id)
except StopIteration:
break
- si on veut récupérer toutes les séquences d'un coup : seqRecordList = list(SeqIO.parse('myFile.fas', 'fasta')) (mais prend de la place en mémoire !)
- on peut lire du fastq : avec SeqIO.parse('myFile.fastq, 'fastq'). La qualité se trouve dans le seqRecord : seqRecord.letter_annotations['phred_quality'] sous forme de liste de score phred (avec les valeurs entières recalées comme il faut en soustrayant 33)
- Pour lire un fichier fastq.gz :
with gzip.open('myFile.fastq.gz', 'rt') as fh:
for seqRecord in SeqIO.parse(fh, 'fastq'):
...
Lecture de séquences indexées :
- récupération de toutes les séquences en mémoire sous forme de dictionnaire dont les clefs sont les ids des séquences : myDict = SeqIO.to_dict(SeqIO.parse('myFile.gb', 'genbank'))
- on peut donner à SeqIO.to_dict un argument key_function qui indique une fonction prenant en argument le SeqRecord et renvoyant la clef à mettre dans le dictionnaire retourné !
- récupération d'un dictionnaire sans stocker toutes les séquences en mémoire : myDict = SeqIO.index('myFile.gb', 'genbank') (ne marche qu'avec un nom de fichier, pas un filehandle) : stocke en mémoire seulement les ids et l'emplacement dans le fichier :
- ça permet de manipuler des grands fichiers (mais par contre, il faut le lire)
- on peut alors accéder au SeqRecord par : myDict[myId]
- on peut aussi accéder au texte brut correspondant : myDict.get_raw(myId)
- pour les fichiers encore plus grands, on peut indexer un ou plusieurs fichiers sur le disque (l'index est une base sqlite3 !) :
- construction de l'index qui est créé sur le disque : ind = SeqIO.index_db('myIndex.idx', 'myFile.gbk', 'genbank') (on peut donner une liste de fichiers à la place d'un seul fichier : ind = SeqIO.index_db('myIndex.idx', ['myFile1.gbk', 'myFile2.gbk'], 'genbank'))
- l'objet ind permet de travailler sur les séquences comme avec un dictionnaire :
- ind.keys() : pour avoir les ids
- ind[myId] : pour accéder à un SeqRecord
- ind.get_raw(myId) : pour accéder au texte brut.
- on peut dans une autre session utiliser directement l'index : ind = SeqIO.index_db('myIndex.idx')
- on peut transformer l'id en passant une fonction qui prend l'id de départ en entrée et renvoie le nouvel id qui doit être unique :
def myCallback(seqName):
return seqName + '_2'
ind = SeqIO.index_db('myIndex.idx', 'myFile.gbk', 'genbank', key_function = myCallback)
Attention, quand on relit l'index, il faut repasser cette même fonction, sinon, ça casse : ind = SeqIO.index_db('myIndex.idx', key_function = myCallback)
- pour les gros fichiers conservés compressés :
- gzip n'est pas adapté à l'accès random, il faut utiliser BGZF (Blocked GNU Zip Format), utiliser bgzip qui vient avec les samtools (ou utiliser le module biopython Bio.bgzf)
- on peut utiliser les fichiers bgz avec SeqIO.index ou SeqIO.index_db exactement pareil : ind = SeqIO.index('myFile.gbk.bgz', 'genbank')
Ecriture de séquences :
- print(seqRecord.format('fasta')) : pour imprimer un SeqRecord au format fasta (idem avec genbank)
- on peut aussi utiliser : SeqIO.write(seqRecordList, 'myFile.fas', 'fasta') (écrase le contenu du fichier s'il existait déjà).
- on peut écrire au format fastq, pourvu que le SeqRecord ait une qualité : seqRecord.format('fastq')
- si on veut écrire dans un fichier : fh.write(seqRecord.format('fasta'))
Conversion d'un fichier fasta et qualité 454 en fichier fastq :
with open(args.outFastq, 'w') as fhOut:
for seqRecord in PairedFastaQualIterator(fhFasta, fhQual):
fhOut.write(seqRecord.format('fastq'))
Récupération de séquences directement sur Entrez :
from Bio import Entrez
from Bio import SeqIO
Entrez.email = "email@example.com"
handle = Entrez.efetch(db="nucleotide", rettype="fasta", retmode="text", id="6273291")
seq_record = SeqIO.read(handle, "fasta")
handle.close()
Récupération de séquences directement sur Uniprot :
from Bio import ExPASy
from Bio import SeqIO
handle = ExPASy.get_sprot_raw("O23729")
seq_record = SeqIO.read(handle, "swiss")
handle.close()
Conversion de format (genbank en fasta par exemple) :
- count = SeqIO.convert('myFile.gbk', 'genbank', 'myFile.fas', 'fasta')
Copyright python-simple.com
programmer en python, tutoriel python, graphes en python, Aymeric Duclert