Lecture et écriture des 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'):
      ...
    

• 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')

• 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'))

    with open(args.outFastq, 'w') as fhOut:
        for seqRecord in PairedFastaQualIterator(fhFasta, fhQual):
            fhOut.write(seqRecord.format('fastq'))
  

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()
  

from Bio import ExPASy
from Bio import SeqIO
handle = ExPASy.get_sprot_raw("O23729")
seq_record = SeqIO.read(handle, "swiss")
handle.close()
  

• count = SeqIO.convert('myFile.gbk', 'genbank', 'myFile.fas', 'fasta')