Manipulation des séquences

Pour savoir quelle version de biopython est utilisée : import Bio; print Bio.__version__ (help(Bio) permet de voir quel est le fichier .py utilisé).

Définition d'une séquence :

• Utilisation : from Bio.Seq import Seq
• mySeq = Seq('AGCGATGTAC') : immutable, comme le type string
• mySeq = mySeq + 'AAAAA' : ajout d'une séquence à la fin
• mySeq = 'AAAAA' + mySeq : ajout d'une séquence au début
• isFound = 'GAT' in mySeq : test si sous-chaîne présente.
• print(mySeq[1]) : origine à 0
• print(len(mySeq)) : longueur de la séquence
• print(str(mySeq)) : conversion en string
• print(mySeq.tostring()) : conversion en string, mais la méthode est découragée (deprecated).
• print(mySeq.reverse_complement()) : reverse complement
• print(mySeq.count('AT')) : compte le nombre de AT (non overlappant)
• isTrue = mySeq.startswith('AGC') || mySeq.endswith('ATA') : regarde si démarre ou fini avec une certaine séquence
• pos = mySeq.find('ATG'); pos = mySeq('ATG', start = 5) : première occurence ou -1 si pas trouvé
• prot = mySeq.translate() : donne un objet Seq protéine (traduction).
• mySeq[2:5] : sous-séquence
• si on définit 2 séquences mySeq1 = Seq('ACGTAC'); mySeq2 = Seq('ACGTAC'), alors mySeq1 == mySeq2 est False, tandis que str(mySeq1) == str(mySeq2) est True
• pour convertir une protéine du code à une lettre en code à 3 lettres : from Bio.SeqUtils import seq3; seq3('ALVYC')

Quand on définit une Seq, il vaut mieux la définir avec un alphabet : from Bio.Alphabet import DNAAlphabet; seq = Seq("ACGATGAC", DNAAlphabet()). C'est nécessaire si par exemple on veut sortir un seqRecord au format genbank. Alphabets possibles :

• from Bio.Alphabet import ProteinAlphabet; seq = Seq('ACGHA', ProteinAlphabet())
• from Bio.Alphabet.IUPAC import IUPACUnambiguousDNA; seq = Seq('ACAGT', IUPACUnambiguousDNA())
• from Bio.Alphabet.IUPAC import IUPACAmbiguousDNA; seq = Seq('AMAHT', IUPACAmbiguousDNA())

Conversion d'une protéine avec code à lettres en code à une lettre :

• faire d'abord import Bio.SeqUtils
• utiliser le dictionnaire :Bio.SeqUtils.IUPACData.protein_letters_3to1

Table de traduction en protéines :

• import Bio.Data.CodonTable
• print(Bio.Data.CodonTable.standard_dna_table) : affiche la table de traduction.
• Bio.Data.CodonTable.standard_dna_table.forward_table : le dictionnaire codon -> acide aminé à une lettre (codons stop absents)

Séquences mutables :

• pour créer une séquence mutable à partir de rien : from Bio.Seq import MutableSeq; mySeq = MutableSeq('ACGTACGT').
• on peut aussi créer une séquence mutable à partir d'une séquence normale qui sera un objet indépendant : mySeq = oldSeq.tomutable().
• mySeq[0:3] : une séquence mutable avec les 3 premières bases, indépendante de la séquence mySeq de départ.
• mySeq[2] = 'A' : change la 3ème base (donne ACATACGT).
• mySeq.reverse_complement() : agit en place (donne ACGTACGT).
• mySeq[0:4] = 'AAAA' : changement de plusieurs bases (donne AAAAACGT).
• mySeq[0:1] = 'AAAAAAA' : on peut remplacer une partie par une partie qui n'a pas la même longueur (donne AAAAAAACGTACGT).
• mySeq[0:2] = '' : détruit les deux premières bases (donne GTACGT).
• mySeq[0:0] = 'CCC' : insère au début (donne CCCACGTACGT).
• mySeq.append('A') : rajoute une base à la fin (on ne peut pas en rajouter plusieurs, donne ici ACGTACGTA).
• mySeq.extend('AAA') : rajoute une séquence à la fin (donne ici ACGTACGTAAA).
• pour créer un seq record avec la qualité : seqRecord = SeqRecord(Seq('GAGCTAT'), id = 'mySeq', letter_annotations = {'phred_quality': [ord(x) - 33 for x in 'GGGG-9-']})
• mySeq.toseq() : renvoie une séquence immutable.

SeqRecord : séquence avec un id, des annotations, des features, ... :

• attributs :
  • import : from Bio.SeqRecord import SeqRecord
  • .seq : l'objet Seq
  • .id : l'id
  • .name : le nom commun, comme par exemple un accession number.
  • .description : la description.
  • .letter_annotations : dictionnaire des annotations par lettre, par exemple scores phred : {'phred_quality': [38, 38, 38, 38, 12, 24, 12]}
  • .annotations : dictionaire des annotations (keywords, comment, taxonomy, ...).
  • .features : la liste des SeqFeature.
  • .dbxrefs : la liste des cross-références aux bases.
• construction : seqRecord = SeqRecord(Seq('ACGATCGATGT'), id = 'mySeq', description = 'test')
• si on ne veut pas mettre de description, utiliser seqRecord = SeqRecord(Seq('ACGATCGATGT'), id = 'mySeq', description = '') (sinon, sur la première ligne du fasta, on aura <unknown description>).
• slicing : print(seqRecord[45:78]) : les features complètement incluses sont conservées seules et décalées comme il faut (attention, les annotations type keywords ne sont pas transférées). Origine à 0 et dernier index exclu, comme habituellement en avec les séquences python.
• le slicing marche aussi avec des séquences type fastq (coupe bien la qualité comme il faut).
• seqRecord.reverse_complement() :
  • renvoie le reverse complement, mais sans modifier l'objet initial.
  • le nouveau SeqRecord renvoyé a un id, un name et une description à unknown qu'il faut alors renseigner (au moins pour l'id !)
• on peut concaténer des seqRecord, les features seront calculées en conséquence : seqRecord = seqRecord1 + seqRecord2 + seqRecord3 (utile par exemple pour concaténer des exons)
• si on veut changer l'id d'un SeqRecord pour l'écrire sur le disque par exemple en format fasta, il faut aussi mettre la description à la chaîne vide : seqRecord.id = 'myNewId'; seqRecord.description = '' (sinon, il va remettre l'ancienne description).

Bio.SeqFeature.SeqFeature :

• attributs :
  • .type : type de feature (par exemple CDS)
  • .location : objet Bio.SeqFeature.FeatureLocation ou Bio.SeqFeature.CompoundLocation
  • .strand : le strand de la feature (1, -1, ou 0 si unknown, ou même None si ça n'a pas de sens)
  • .qualifiers : dictionnaire des qualifiers
• myFeature.extract(mySeq) : renvoie une extraction de mySeq correspondant à la feature (si mySeq est un SeqRecord, renvoie un SeqRecord, si mySeq est un Seq, renvoie un Seq)
• on va faire souvent : mySeqRecord.extract(mySeqRecord.features[3])
• myFeature.translate(mySeq) : renvoie la traduction de mySeq correspondant à la feature (si mySeq est un SeqRecord, renvoie un SeqRecord, si mySeq est un Seq, renvoie un Seq)

FeatureLocation :

• import : from Bio.SeqFeature import FeatureLocation
• attributs :
  • .start, .end : les limites de la feature : attention, origine = 0, et end est la position après la dernière (comme d'habitude en python)> Donc une feature qui est notée 1..10 dans le fichier d'origine va avoir start = 0 et end = 10 !
  • .strand : 1, -1 ou 0
• fonctions :
  • str(myLoc) : représentation textuelle
  • myLoc.extract(mySeq) : extrait la séquence correspondant à la feature sur la séquence qui porte cette feature.
  • myLoc.extract(mySeqRecord) : marche aussi avec un SeqRecord et renvoie alors un seqRecord.
• il y a aussi un attribut .parts qui renvoie une liste avec un seul élément, lui-même (pour compatibilité avec CompoundLocation).

CompoundLocation :

• représente une localisation avec plusieurs segments (par exemple, la location d'un ARNm)
• import : from Bio.SeqFeature import CompoundLocation
• a les mêmes attributs que FeatureLocation : .start, .end, .strand et la même fonction extract()
• mais attention : ce n'est pas une classe dérivée de FeatureLocation.
• pour avoir les segments de la location, faire mylocation.parts, cela renvoie une liste de FeatureLocation.

Reconstruire un mRNA à partir des features de type exon d'un SeqRecord :

exonFeatures = filter(lambda x: x.type == 'exon', seqRecord.features)
newSeqRecord = reduce(lambda x, y: x + y, [feat.location.extract(seqRecord) for feat in exonFeatures])
  

conversion code 1 lettre / code 3 lettres protéique :

• from Bio.SeqUtils import seq3; print(seq3('MAI'))
• from Bio.SeqUtils import seq1; print(seq1('MetAlaIle'))