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