Licence CC BY-NC-ND François Rechenmann & Thierry Parmentelat

Recalage par recherche des RBS

Dans ce complément, nous allons implémenter l'algorithme de recalage par RBS.

Variante des régions codantes, avec RBS

Afin d'implémenter l'idée de recalage telle qu'elle est décrite dans la vidéo, nous allons modifier coding_regions_one_phase comme suit. Vous remarquerez que la différence par rapport à la version simple est très faible :

In [1]:
# on importe les fonctions annexes que nous avions utilisées
# pour écrire la version simple de cet algorithme
from w3_s03_c2_next_codon import next_start_codon, next_stop_codon

# avec à nouveau une longueur minimale de 300 par défaut
# et en utilisant le recalage RBS
def coding_regions_one_phase_rbs(dna, phase, rbs, minimal_length=300):
    # on initialise index à la phase; avec next_start_codon
    # et next_stop_codon on reste toujours sur la même phase
    index = phase
    # les résultats sont retournés sous la forme d'une liste 
    # de couples [start_gene, stop_gene]
    genes = []
    # on calcule stop1 qui désigne pour nous le stop "de gauche"
    # à ce stade, c'est le premier STOP à partir de la phase
    stop1 = next_stop_codon(dna, index)
    # s'il n'y a pas du tout de stop dans toute la séquence
    # on a terminé
    if not stop1:
        return genes
    # la boucle principale
    while True:
        # on cherche le premier STOP à partir de stop1
        # pour trouver le stop "de droite"
        stop2 = next_stop_codon(dna, stop1+3)
        # s'il n'y a plus de stop, on a fini
        if not stop2:
            return genes
        # mais il faut qu'il soit assez loin
        if stop2 - stop1 < minimal_length:
            # c'est trop court, on saute ce fragment
            stop1 = stop2
            continue
        # à ce stade on a trouvé un ORF, reste à trouver le bon start
        start = next_start_codon(dna, stop1)
        # s'il n'y en a pas: c'est qu'on ne trouvera plus rien
        # et donc on a fini
        if not start:
            return genes
        # regardons si on peut trouver un RBS dans la region codante
        next_rbs = dna.find(rbs, start)
        # est-ce bien dans la region ?
        if start <= next_rbs <= stop2:
            # le RBS est bien dans la région:
            # on recale le debut de la région codante comme 
            # le prochain START à droite du RBS
            start = next_start_codon(dna, next_rbs)
        if start and stop2 - start < minimal_length:
            # si la region est trop petite, on l'ignore
            pass
        else:
            # cette fois on a trouvé un gène, on l'ajoute dans les résultats
            genes.append( [start, stop2] )
        # on peut passer à l'ORF suivant
        stop1 = stop2

the rest of 25 divided by 10 is 5
dealing with 1
dealing with 2
dealing with 4
PHASE 0
found at index 0 ATG
found at index 9 ATG
found at index 54 ATG
PHASE 1
found at index 40 ATG
PHASE 2
found at index 5 ATG
found at index 32 ATG
ATGCGATGTATGCGTGCAGCTGCTAGCTCGTAATGTCGTCATGGATCATCGATCATGG
PHASE 0
found at index 30 TAA
PHASE 1
PHASE 2
found at index 23 TAG

Des données réelles

Nous allons voir le résultat des deux versions de regions_codantes sur une bactérie réelle, la célèbre Escherichia Coli, que vous pouvez consulter sur ENA sous la clé U00096, mais que nous avons téléchargé pour vous :

In [2]:
with open("data/escherichia-coli-U00096") as input:
    e_coli = input.read()
print("Escheria Coli a {} nucléotides".format(len(e_coli)))

Escheria Coli a 4641652 nucléotides

Si on importe l'algorithme simple :

In [3]:
# la recherche de régions codantes sur une phase
# telle que nous l'avons vue dans la séquence 2
from w3_s02_c1_coding_regions_v1 import coding_regions_one_phase

counter =  2
counter =  4
counter =  8
counter =  16
counter =  32
counter =  64
after the loop
subtilis has 4215367 bases
We found 788 genes on phase 0
average gene length 830.2119289340102
min = 300, max = 7629
Percentage of coding region 0.1551957397778177

Il nous donne avec cet exemple :

In [4]:
for phase in 0, 1, 2:
    algo_simple = coding_regions_one_phase(e_coli, phase)
    # calculons la taille moyenne
    how_many = len(algo_simple)
    total_length = sum ( stop-start for start, stop in algo_simple )
    average_size = total_length / how_many
    print("L'algorithme simple trouve {} regions, taille moyenne = {:.02f} sur la phase {}"
          .format(how_many, average_size, phase))

L'algorithme simple trouve 983 regions, taille moyenne = 798.41 sur la phase 0
L'algorithme simple trouve 1016 regions, taille moyenne = 822.90 sur la phase 1
L'algorithme simple trouve 987 regions, taille moyenne = 819.36 sur la phase 2

Si maintenant nous utilisons la séquence AGGAGG comme RBS, nous obtenons cette fois :

In [6]:
rbs_coli = 'AGGAGG'
for phase in 0, 1, 2:
    algo_rbs = coding_regions_one_phase_rbs(e_coli, phase, rbs_coli)
    # calculons la taille moyenne
    how_many = len(algo_rbs)
    total_length = sum ( stop-start for start, stop in algo_rbs )
    average_size = total_length / how_many
    print("L'algorithme RBS trouve {} regions, taille moyenne = {:.02f} sur la phase {}"
          .format(how_many, average_size, phase))

L'algorithme RBS trouve 965 regions, taille moyenne = 791.64 sur la phase 0
L'algorithme RBS trouve 990 regions, taille moyenne = 811.70 sur la phase 1
L'algorithme RBS trouve 964 regions, taille moyenne = 802.97 sur la phase 2
In [7]:
from w3_s03_c2_next_codon import next_start_codon, next_stop_codon

Vous pouvez donc voir qu'on obtient avec cette méthode un recalage assez marginal, qui conserve les ordres de grandeur, puisqu'on observe :

  • sensiblement le même nombre de régions codantes ; la différence s'explique par des régions où le recalage raccourcit la séquence au dessous du seuil de longueur.
  • des tailles de régions codantes là encore dans le même ordre de grandeur.