Couper une séquence nucléotidique en deux

Version imprimable

Bonjour,
Alors voilà comme dit dans le titre mon objectif est de couper une séquence nucléotidique (chaîne de caractères avec que des A,T,G ou C). Le problème c'est qu'il faudrait couper en deux la séquence mais par un multiple de 3 afin de ne pas couper un codon (représenté par 3 de ces 4 lettres). De plus, j'aimerais ne pas couper des suites de 'AAA' ou 'AAG' ou les deux en même temps, s'ils se retrouvent au milieu de la séquence. En revanche, s'il n'y a pas ces lettres au milieu de la séquence, il faudrait la couper en deux.
Voici une partie de mon code qui sert à regarder si un AAA, AAG ou les 2 en même temps sont au milieu de la séquence :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 if 'AAA' in seq[i-2:i+3] or 'AAG' in seq[i-2:i+3]: l = 0 for i in range(i-2, i+1): j = i while True: if seq[i:j+3] in s: j += 3 else: if j > l: l = j break print(str(seq[:l])) else: if len(seq) % 2 == 0: print(seq[:len(seq) // 2]) else: print(seq[:(len(seq) + 4) // 2])
Input :
ATGTGAGTCAAAAAATGCCTATGA
Output :
ATGTGAGTC AAAAAATGCCTATGA ou ATGTGAGTCAAAAAA TGCCTATGA
La longueur des suites successives de AAA ou AAG varie en fonction des séquences, on peut trouver des AAAAAA, des AAGAAG ou des AAAAAG (ou AAGAAA).

Salut,

Un truc comme çà peut être:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
>>> seq = 'ATGTGAGTCAAAAAATGCCTATGA'
>>> frag = []
>>> cutoff = False
>>> for i in range(0, len(seq), 3):
...     s = seq[i:i+3]
...     if s in ('AAA', 'AAG'):
...         frag.append(s)
...         cutoff = True
...     else:
...         if cutoff:
...             print (''.join(frag))
...             cutoff = False
...             frag = []
...         frag.append(s)
...
ATGTGAGTCAAAAAA
>>> print(''.join(frag))
TGCCTATGA
>>>

- W

24/09/2018, 21h36
glageosa

Merci beaucoup, je crois que cela devrait le faire.
24/09/2018, 21h59
glageosa

Ah non en fait c'est pas tout à fait ça excuse moi, l'idée c'est de juste couper la séquence en deux et s'il y a un ployA alors il faut le mettre de coté, là ton algo coupe à chaque fois qu'il en voit un.

salut,

un truc dans ce style fait peut-être l'affaire :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
import re
 
def cut (seq, liste):
   middle = int(len(seq)/6)*3                            # le milieu de la sequence, aligné sur 3 nucléotides
   reg = re.compile(r'({})+'.format(r'|'.join(liste)))   # regex créée dynamiquement pour capter les codons à ne pas séparer
   for i in reg.finditer(seq):
      start = i.start()
      end = i.end()
      if start < middle < end:                           # le milieu calculé tombe en plein dans une suite de codons
         if (middle - start) <= (end - middle):          # on équilibre la taille des sous-séquences à renvoyer
            return (seq[:start], seq[start:])
         else:
            return (s[:end], s[end:])
   return (seq[:middle], seq[middle:])

on détermine le milieu (aligné sur 3 nucléotides) de la séquence initiale, et on compare ensuite avec la liste des (début,fin) de chaque suite de codons à ne pas séparer dans la séquence ("AAAAATAAA" par exemple)
si le milieu de la séquence qu'on a calculé tombe au milieu d'une séquence de codons, alors on s'arrange pour équilibrer le plus possible la taille des sous-séquences retournées

et le résultat :

Code:

1
2
3
4
5
6
7
8
>>> liste = ['AAA', 'AAT']
>>> cut('', liste)                              -->      ('', '')
>>> cut('AAA', liste)                           -->      ('', 'AAA')
>>> cut('TTC', liste)                           -->      ('', 'TTC')
>>> cut('TTCGAC', liste)                        -->      ('TTC', 'GAC')
>>> cut('AATAAA', liste)                        -->      ('', 'AATAAA')
>>> cut('ATGTGAGTCAAAAAATGCCTATGA', liste)      -->      ('ATGTGAGTC', 'AAAAAATGCCTATGA')
>>> cut('ATGTGAAAAAAAGTCTGCCTATGA', liste)      -->      ('ATGTGAAAAAAA', 'GTCTGCCTATGA')

25/09/2018, 19h27
glageosa

Merci tu as très bien compris ce que je voulais faire, cependant il reste un minuscule hic, est-ce qu'il cherche la regex en prenant compte du fait que c'est 3 par 3 ? Parce que j'ai un problème avec une séquence, il coupe 2 fois, une fois mal en voyant un AAG mais qui n'est pas dans le cadre de lecture (c'est-à-dire pas multiple de 3) et une dernière fois bien. Je m'en suis rendu compte en mettant un print à la place du return, et lorsque l'on met return il ne coupe cette séquence qu'une seule fois mais de la mauvaise façon. J'ai ensuite regardé l'exécution sur python tutor parce que je ne comprenais pas et j'ai remarqué que la première fois lorsqu'il coupe mal, il execute l'étape de rééquilibrage (print(seq[:start], seq[start:])) et quand il coupe bien il execute le dernier print (print(seq[:middle], seq[middle:]))
Merci à toi.
26/09/2018, 21h35
BufferBob

Citation:

Envoyé par glageosa

est-ce qu'il cherche la regex en prenant compte du fait que c'est 3 par 3 ?

non t'as complètement raison c'est ma regex qui ne fait pas le job comme prévu :oops:
j'aurai aimé trouver une solution reposant le plus possible sur le moteur d'expressions rationnelles (déjà existant donc, plutôt que d'avoir à réimplémenter en tout ou partie un algo de recherche), mais rien, peau de balle.
27/09/2018, 10h50
disedorgue
Bonjour,

Est-ce que ce ne serait pas plus simple de transformer d'abord ce que tu ne dois pas couper, par exemple:
- AAAAAA devient I
- AAGAAA devient J
- AAAAAG devient K
- AAGAAA devient L
- AAA devient M
- AAG devient N
Ensuite tu coupes et tu refait la transformée inverse ?

PS: Un petit détail, sur la transformée initiale, il faut la faire en traitant les grandes occurrences d'abord et une fois fait, on traite les petite (de 3 lettres).

Bon, après, j'ai pas tout compris, comme par exemple, peut-on couper une séquence du genre 'AAGAAGAAAAAA' en 'AAGAAG' et 'AAAAAA' ?
Si ce n'est pas le cas, ce que je propose ne fonctionnera pas.

En s'apuyant sur les positions des codons à ne pas couper, ça ne peut pas fonctionner ?

Du genre :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
import re
 
def scinder(seq) :
    ex = list()
    lseq = len(seq)
    for s in re.finditer("AA(?:A|G)", seq) :
        if not s.start() % 3 :
            ex.append(s.start())
    mm = lseq // 2 - lseq % 2
    for i in range(mm, lseq, 3) :
        if i not in ex :
            break
    return (seq[0:i], seq[i:])
 
seqs = (
    'ATGTGAGTCAAAAAATGCCTATGA',
    'ATGTGAGTCAAGAAATGCCTATGA',
    'ATGTGAGTCAAGAAATGCCTATGA',
    'ATGTGAGTCATGTGTTGCCTATGA',
    'ATGTGAAAAAAAGTCTGCCTATGA',
)
 
for s in seqs :
    parts = scinder(s)
    print(s, parts, len(parts[0]), len(parts[1]))
    # Partie insécable dans la seconde moitié ?
    m = re.search("(?:AA(?:A|G))+$", parts[0])
    if m :
        ip = (parts[0][0:m.start()], parts[0][m.start():] + parts[1])
        print(s, ip, len(ip[0]), len(ip[1]))
    print('\n')

Edit : Oui en effet, cela ne marche pas bien, alors qu'au final avec un truc plus simple, ça fonctionne.

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
def scinder(seq, insecables) :
    lseq = len(seq)
    # Dict pas vraiment nécessaire, mais permet plus de clarté
    dd = {i:seq[i:i+3] in insecables for i in range(0, lseq, 3)}
    mm = lseq // 2 - lseq % 2
    for i in range(mm, lseq, 3) :
        if not dd[i] :
            break
    else :
        # Impossible de scinder la séquence entre mm et fin de séquence
        # Alors tentative de la scinder entre mm et début de séquence
        for i in range(mm, 0, -3) :
            if not dd[i] :
                i += 3
                break
        else :
            # La séquence entière est insécable
            i = lseq
 
    return (seq[0:i], seq[i:])

en revoyant ma copie j'ai ça :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 def cut (seq, liste): middle = int(len(seq)/6)*3 # le milieu de la séquence, alignée sur 3 nucléotides s1, s2 = seq[:middle], seq[middle:] # deux parties s1 et s2 à peu près égales r1 = re.search(r'((?:' + r'|'.join(liste) + r')+)$', s1) # match les codons insécables en fin de s1 r2 = re.search(r'^((?:' + r'|'.join(liste) + r')+)', s2) # match les codons insécables en debut de s2 if (r1 is not None) and (r2 is not None): # si les deux match, c'est qu'on est au milieu d'une séquence insécable if len(r1.group()) < len(r2.group()): # selon de quel coté est la partie insécable la plus courte... i, _ = r1.span() s1, s2 = s1[:i], s1[i:] + s2 # ...on reaffecte s1 et s2 else: _, i = r2.span() s1, s2 = s1 + s2[:i], s2[i:] # ...on reaffecte s1 et s2 return (s1, s2)
normalement ça devrait fonctionner comme attendu, à valider néanmoins :aie:

le principe est finalement très simple; on coupe la séquence initiale en deux parties s1 et s2 à peu près égales, puis on cherche les suites de codons à la fin de s1 et au début de s2
si l'une des deux regex ne remonte rien, alors on peut renvoyer s1 et s2 telles quelles, sinon on cherche à équilibrer leurs tailles en faisant passer la plus petite suite insécable trouvée dans l'autre chaine

j'ai pas testé les approches des petits copains mais c'est probablement tout aussi valable (soit se simplifier la vie en se ramenant à des suites de caractères plutôt que de jongler avec des multiples de 3, soit à l'inverse faire de l'arithmétique de base sur la séquence pour retrouver ses petits)