Discussion :

[Totoprosper]

Bonjour,

Je suis un vrai débutant en matière de programmation. J’utilise les langages de programmation awk, grep, sed et python pour analyser des séquences d’ADN. Le problème auquel je fais fasse depuis un certain moment est le suivant :
J’ai 1024 fichiers en format fasta. Ils sont du genre :
…….BamTo_AAAAA.fasta
…….BamTo_CCCCC.fasta
…….BamTo_GGGGG.fasta
…….BamTo_TTTTT.fasta
…….BamTo_ACAAA.fasta
…….BamTo_AGAAA.fasta
…….BamTo_ATAAA.fasta
……………………….. Avec A,C,G,T les 4 constituants (bases) des séquences d’ADN (ici 45 = 1024 = nombre de fichiers)

Chaque fichier contient à peu près 200 à 300 séquences avec des tailles différentes pouvant aller de 50 à 240 constituants (bases).

Dans chaque fichier, il y a des séquences répétées 1x, 2x, 2x, 3x, 4x, 5x,…..Nx. Le fichier …….BamTo_AAAAA.fasta peut par exemple contenir les séquences suivantes:

----------------------------------------------G----- séquence répétée 1x

------------------------------T--------- séquence répétée 2x
------------------------------T---------

------------------------------------A--------- séquence répétée 2x (avec une erreur)
------------------------------------G---------

---------------------------------C---------------------------- séquence répétée 3x (sans erreur)
---------------------------------C----------------------------
---------------------------------C----------------------------

-----------------------------------------C--------------------------- séquence répétée 3x (avec une erreur)
-----------------------------------------A---------------------------
-----------------------------------------A---------------------------

------G-----------------------------------A-------------------------- (séquences répétée 5x avec 3 erreurs)
------G-----------------------------------A--------------------------
------C-----------------------------------G--------------------------
------T-----------------------------------A--------------------------
------G----------------------------------A---------------------------

J’aimerais maintenant créer un fichier consensus (consensus_fichier1.fasta) qui contient : toutes les séquences répétées 1x et 2x de tous les fichiers, un représentant de chaque séquence répétée au moins 3x (après correction des erreurs). Dans le fichier consensus, une même séquence peut être répétée xfois, cela ne pose aucun problème.

Egalement, je souhaiterais un deuxième fichier consensus (consensus_fichier2.fasta) similaire au premier mais qui se présente comme suit :

Sequence1_AAAAA --------------------------------T-----------------------
Seuence2_ACAAA ---------------------------------------------------G----------------
Sequence3_AACAA -------------------------------A----------------

Chaque ligne du deuxième fichier consensus contient le nom de la séquence (sequence n_NNNNN ; {NNNNN} étant le suffix du fichier d’origine d’où provient la séquence n) suivi de tabulation et la séquence elle-même.

Avec le plus grand respect, quelqu’un pourrait-il m’aider avec quelques codes de awk, sed, grep ou de python à accomplir ce travail ?

D’avance, mille mercis

[N_BaH]

je vais faire ma mauvaise tête, mais "répéter 1 fois", tel que tu le montres, ce n'est pas possible.
répéter, c'est re-dire.
Or, ----------------------------------------------G-----, est présent une seule fois, donc n'est pas répété !
par contre, ------------------------------T---------, est effectivement répété une fois : la deuxième ligne répète la première.

je suis perplexe devant le sujet : autant je peux identifier des lignes identiques, autant comment puis-je ne pas considérer qu'une ligne, qui ressemble, qui comporte (selon tes termes) une/des erreur(s), n'est pas une nouvelle séquence ?

il me faudrait quelques extraits de fichiers représentatifs, et la sortie attendue.

[Flodelarab]

Bonjour

L'ordre des lignes a-t-elle une importance ?
La phrase à droite de la séquence est-elle dans le fichier ou seulement sur le forum ?
Les tirets signifient-ils "n'importe quel caractère" ou y a-t-il vraiment des tirets dans le fichier ?
Et même questions que N_BaH.

[Totoprosper]

Bonjour,
Merci N_BAH pour votre intervention. J'aimerais tout d'abord rectifier que s'il y a dans un fichier .....BamTo_NNNNN.fasta une séquence répétée 2x avec une identité de 100%, seule une des deux séquences est présente dans le fichier consensus. S'il y a une erreur dans l'une des séquences, les deux sont prises dans le fichier consensus.

Les erreurs proviennent du séquençage et de l'amplification des séquences. S'il y a une similitude d'au moins 95%, entre deux séquences, elles peuvent être
considérées mêmes séquences, la présence d'une troisième séquence va permettre de corriger justement les erreurs.

Jai joint deux fichiers type et le résultat attendu (fichier consensus). J'ai mis des lettres et des numéros pour des raisons de clarification, mais les fichiers originaux n'en disposent pas. Merci

[Totoprosper]

Bonjour,

Merci Flodelarab. L'ordre des lignes n'a aucune importance. Le nom des séquences à droite ne se trouve pas dans les fichiers. C'était pour
clarifier ma question. Il n'y a pas non plus de tirets dans les fichiers (les tirets sont dans les nom des fichiers.). Dans les fichiers, chaque ligne est une séquence. Pas de ligne entre séquences.

Merci beaucoup

[Flodelarab]

La question laisse perplexe. A priori, on fait ceci :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
32
33
34
35
36
37
38
39
40
41
42
43
$ ls adn_[GT]*.txt
adn_GGAAT.txt  adn_TTTCC.txt
$ sort adn_[GT]*.txt |uniq -c
      2 AACCTCAATGCTCTGCTACCCTGAATTCTCAGGTTGAGCTGGAGGTTCCCCAGCTCTGCAGCTTCATCCTGAAGACCAGCCAGTGCACCCTGAAAGAGGTGTATGGCTTCAACCCTGAGGGCAAAGCCTTGCTGAAGAAAACCAAGAACTCTGAGGAGTTTG
      3 AGGTGCCCGCCACCACGCCCGGAACTTTTTGTATTTTCAGTAGAGACGGGGTTTCACCTTGGTCTCGATCTCCTGACCTCGTGATCCACCCGCCTCAGCCTCCCAAAGTGCTGGGATTACAGGCGTGAGCCACTGCGCCCGGCCTCACC
      2 AGGTGCCCGCCACCACGCCCGGCCAATTTTTGTATTTTTAGTAGAAGACATGGTTTCACCATGTTGGCCAGGCTGATCTCGATCTCCTAACCTCAAGTGATCCACCCACCTCGGCCTCCCAAAGTGCTGGGTTTACAGGTGTGAGCCACTGCGCCCGGCCTCACC
      1 AGGTGCCCGCCACCACGCCCGGCTAATTTTGTATTTTTAGTAGAGACGGGGTTTCTCCATGTTGGTCAGGGTGGTCTCGAACTCCCAACCTCAGGTGATCCGCCCGCCTCCGCCTCCCAAAGTGTTGGAATTACAGATGTGAGCCACTGCGCCCGGCCTCACC
      1 AGGTGCCCGCCACCACGCCCGGCTAATTTTTCTGTATTTTTGTAGAGACAGGGTTTCACTGTGTTAGCCAGGATGGTCTCGATCTCCTGACCTCGTGATCTGCCTGCCTCGGCCTCCCAAAGTGCTGGGATTACAGGCATGAGCCACTGCGCCCGGCCTCACC
      1 AGGTGCCCGCCACCACGCCCGGCTAATTTTTGTATTTTTAGTAGAGACAGGGTTTCATCATGTTGGCCAGGTGGGTCTCAAACTCCTGACCTCAGGTGGTCCTTCTGCCTTGGCCTCCCAAAGTGCTGGGATTATAGGCGTGAGCCACTGCGCCCGGCCTCACC
      1 AGGTGCCCGCCACCACGCCCGGCTAATTTTTTGTATTTTTAGCAGAGTCGGGGTTTCACCATGTTAGCCAGGATGGTCTCGATCTCCTGACCTCCTAATCCGCCCGCCTCAGCCTCCCAAAGTGCTGGGATTACAGGCATGAGCCACTGCGCCCGGCCTCACC
      1 AGGTGCCCGCCACCACGCCCGGCTAATTTTTTGTATTTTTAGTAGAGACGGGATTTCACCGTGTTAGCCAGGATGGTCTCGATCTCCTGACCTGGTGATCCGCCCGCCTCGGCCTCCCAAAGTGCTAGGATTATAGGTGTGAACCACTGCGCCCGGCCTCACC
      1 AGGTGCCCGCCACCACGCCCGGCTAATTTTTTGTATTTTTAGTAGAGACGGGATTTCACCGTGTTAGCCAGGATGGTCTCGATCTCCTGACCTTGTGATCCGCCCGCCGCGGACTCCCAGAGTGCTGGGATTGCAGGCGAGAGCCACTGCGCCCGGCCTCACC
      1 AGGTGCCCGCCACCACGCCCGGCTAATTTTTTGTATTTTTAGTAGAGACGGGGTTTCACCATGTTAGCCAGGATGGTCTCGATATCCTGACCTTGTGATCTGCCCGCCTTGGCCTCCCAAAGTGCTGGGATTACAGGCTTGAGCCACTGCGCCCGGCCTCACC
      1 AGGTGCCCGCCACCACGCCCGGCTAATTTTTTGTATTTTTTGTAGAGATAGGATTTCACTGTGTTAGCCAGGATGGTCTCGATCTCCTGACCTCGTTATCCACCCACCTGAACCTCCCAAAGTGCTGGGATTACAGGCTTGAGCCACTGCGCCCGGCCTCACC
      1 AGGTGCCCGCCACCACGCCCGGCTAATTTTTTTGAATTTTTAGTACAGACGGGGTTTACCGTGTTAGCCAGGATGGTCCCCTTCTCCTGACCTCGTGATCCACCCGCCTCGGCCTCCCAAAGTGCTGGGATTACAGGCGTGAGCCACTGCGCCCGGCCTCACC
      1 AGGTGCCCGCCACCACGCCCGGCTAATTTTTTTTGTATTTTTAGTAGAGACGGGGTTTCACCACGGTCTCGATCTCCTGACCTCGTGATCCACCCGCCTTGGCCTCTCAAAGTGCTGGGATTAAAGGTGTGAGCCACTGCGCCCGGCCTCACC
      1 AGGTGCCCGCCACCACGCCCGGCTATTTTTTGTATTTTTAGTAGAGACGGGGTTTCACCGTGTTAGCCAGGATGGTCTCGATCTCCTGACCTCCTGATTCACCCGTCTCGGCCTCCCGAGGTGCTGGATTTACAGGAGTGAGCCACTGCGCCCGCCTCACC
      1 AGGTGCCCGCCACCACGTTAGGCTAAGTTTTTGTATTTTTAGTAGAGACGGGGTTTCACCACATTAGCCAGGATGGTCTTGATCTCCTGACCTCGTGATCCACCTGCCTCGGCCTCCTAAAATGCTGGGATTACAGGCGTGAGCCACTGCGCCCGGCCTCACC
      1 CACTCCAGAGAGCCCACTACCGTCGGATCGTGCGTGTATCCTACGATCTCGTATGCCGTCTTCTGCTTGAAAAAAAAAAAAAACAAAAATAAAAAAAAAAAATAAAAAAAAAAAAAAAATAAAAAAAAAAATATATTTTTTTTTTTTAATGATACGGCGACCACCGAGATCTACACATACGAGATCCGTAATCGGGAAGCTGAAGCACTCCAG
      1 CAGAAGCCAGAAGTGAGATCCTCGCCCACTGGTCGCCTGCCAAACTGCTTCTCCAAATGGACTCATCTGCTACAGCTTATGGCTCCACAGTTTCCAAGAGGGTGGCATGGCATTATGGTATGTGTCTCTTCCCCTGTGTGAGCACTTCCAAAGTAATGCAGG
      3 CATTATGCACTGTCTTCATTGTAATGTCTTTAGATTAGAGCTGGAACTGAGGCCTACAGGAGAGATTGAGCAGTATTCTGTCAGCGCAACCTATGAGCTCCAGAGAGAGGACAGAGCCTTGGTGGATACCCTGAAGTTTGTAACTCAAGCAGAAGGTGAGTA
      1 CCAAGGATGGTCTTTGATTTCAGGAACAGTGAACCCTTGCTCTACCAATGCTTTCATACGTTTAGCCCAATCTTGGATAGAATATTGCTCTGCAAAGTCAGTAAGGTTCTTAGCAGCAAGAGTCCACCAATCAGAAATGTAGGTGACAAGTGCGCTATAAACC
      1 CCAAGGATGGTCTTTGATTTCAGGAACAGTGAACCCTTGCTCTACCAATGCTTTCATACGTTTAGCCCAATCTTGGATAGAATATTGCTCTGCAAAGTCAGTAAGGTTCTTAGCAGCAAGAGTCCACCAATCAGAAATGTAGGTGACAAGTGTGCTATAAACC
      1 CCAGGATGGTCTCCGATCTCCTGACCTCGTGATCCGCCTGCCTTGGCCTCCCGAAGTGCTGGGATTACAGGCGTGAGCCACTGCGCCCGGCCTCACC
      1 CCAGGATGGTCTCGATCTCCTGACCTCGTGATCCGCCTGCCTCGGCCTCCGAAAGTGCTGAGATCACAGGCATGAGCCACTGCGCCCGGCCTCACC
      4 CCAGGATGGTCTCGATCTCCTGACCTCGTGATCCGCCTGCCTCGGCCTCCGAAAGTGCTGAGATTACAGGCATGAGCCACTGCGCCCGGCCTCACC
      1 CCTGATGCCATTTTCTTGTCTGATTTTCTTTCTCAGATTGGCTTGGAAGGAAAAGGCTTTGAGCCAACATTGGAAGCTCTTTTTGGGAAGCAAGGATTTTTCCCAGACAGTGTCAACAAAGCTTTGTACTGGGTTAATGGTCAAGTTCCTGATGGTGTCTCT
      7 CCTGTAGAGTGGGAGGCCCTCAGGGACCCGGGTGTAGGAGAGTGCACGGGGCTGGGCGCCCTTCCACGCCCCATGCGCAGATGCCTTACTTGGACAGACCAGGCTGACATTTTCCAGCATTTCCTCTTCTGTAAGACAGGAGAAAGAAATCTGTGAGCTTCC
      1 CTCCAACTCCTGAACTCAAGTGATCGACCCACCTCGGCCTCCCAAAGTGCTGGGATTACAGGTGTGAGCCACTGCGCCCGGCCTCACC
      1 GAACAGGGCTGGGGGAAAGCTGTGGGCTCTAGGTCCCTCCTGCCTGCATCCTCCATACCTTGCAGTTGATCCTGGTGGCACTTCTTGAATCAGCAGTCCCAGGGACTCCACTGGAACTCTCAGCCTCATAGTTGTATGTGTACTTCCGGAGGTGCTTGAATC
      1 GCTAACAAGGTGAAACCCCGTCTCTACTAAAAATACAAAAAATTAGCCGGGCGCGGTGGCGGGCGCCTGTAGTCCCAGCTACTCGGGAGGCTGAGGCAGGAGAATGGTGTGAAC
      2 GGAGGCCGAGGCGGGTGGATCACAAGGTCAGGAGATAGAGACCATCCTGGCTAACACGGTGAAACCCCATCTCTGCTAAAAATACAAAAAATTAGCCAGGCTTGGTGGTGGAAGCCTGTAGTCCCAGCTACTTGGGAGGCTGAGGCAGGAGAATGGTGTGAAC
      1 GGAGGCCGAGGCGGGTGGATCATGAGGTCAGGAGATCAAGACCATCCTGGCTAACACGGTGAAACCCCGTCTCTACTAAAAATATAAAAAATTAGCCAGGCGTCGTGGTGGGCGCCTGTAGTTCCAGCTACTCGTGAGGCTGAGGCAGGAGAATGGTGTGAAC
      1 GGAGGCCGAGGCGGGTGGATCATGAGGTCAGGAGATCAAGACCATCCTGGCTAACATGGTGCATCAGCAGCCC
      1 GGAGGCCGAGGCGGGTGGATCATGAGGTCAGGAGATCGAGACCATCCTGGCTAACAAGGTGAAACCCCGTCTCTACTAAAAATACAAAAAATTAGCCGGGCGCGGTGGCGGGCGCCTGTAGTCCCAGCTACTCGGGAGGCTGAGGCAGGAGAATGGTGTGAAC
      1 GGAGGCCGAGGCGGGTGGATCATGAGGTCAGGAGATCGAGACCATCCTGGCTAACAAGGTGAAACCCCGTCTCTACTAAAAATACAAAAAATTAGCCGGGTGTGGTGGCGGGTGCCTGTAGTCCCAGCTATTTGGGAGGCTGAGGCGGGAGAATGGTGTGAAC
      1 GGAGGCCGAGGCGGGTGGATCATGAGGTCAGGAGATTGAGACCATCCTGGCTAACATGGTAAAACCCCATCTCTACTAAAAATACAAAAAATTTAGCCAGGCATGGTGGCGGGCGCCTGTAGTCCCAGCTACTTGGGAGGCTGAGGCAGGAGAATGGTGTGAAC
      1 GTGAGAGCAGTCAGTTTCTCCTTGGCATGTGAAACTTGTCTCTCCCAATTGAATGAATTCAGATAATCATTAGCTTGCTGTGGGAGTTTTCCCAGGGCTGCTCTGTATTTTCTTACAAATTGATCAATATTGATGTGCTTCAGGTTTCTCTGTACGTTTTCC
      1 TCAAACACAGTGGCAAGTTTACACACAGAAAAAAATACACTGGAGCTTAGTAATGGAGTGATTGTCAAGATAAACAATCAGCTTACCCTGGATAGCAACACTAAATACTTCCACAAATTGAACATCCCCAAACTGGACTTCTCTAGTCAGGCTGACCTGCGC
      2 TCCTCTCAAAATATTCTTGCAAGGTCTCAAATGACATGATGGGCTCACTGA
      1 TCTCTGGTAACTACTTTAAACAGTGACCTGAAATACAATGCTCTGGATCTCACCAACAATGGGAAACTACGGCTAGAACCCCTGAAGCTGCATGTGGCTGGTAACCTAAAAGGAGCCTACCAAAATAATGAAATAAAACACATCTATGCCATCTCTTCTGCT
      1 TGTGTTCCCTCATCTGAGAATCTGGGGCAGGCCCATTTCCATGACCCTTTTCCAGAAGAAGTCCATGCTATGTGGCCAGCTTTCAACAGTGTCTTGATCTCGTTGCGCAGGTCAGCCTGACTAGAGAAGTCCAGTTTGGGGATGTTCAATTTGTGGAAGTAT
      1 TTCTAGGTTTCCTCCAGTGGTGGGGAGGCGGACGAGGAAAAGCAGAGTAGTGGACAGGACACGTCATGTTCCTCATAGACATCTGGAACCTTCTCATGCATCGTTTCCTTCCCTTCTAGGCACTGCCCATTTCCCACCGGAAACTTCAGCCAGCGCTCGCTG

À gauche, la quantité. À droite la séquence nécessairement unique.

Il faut trier avec "sort" car "uniq" le demande.

Était-ce la question ?

[papajoker]

d'après ce que j'ai compris ... il doit être possible d'avoir un nombre X de caractères différents ("95% uniquement doit être unique")
Donc c'est plutôt une question algorithmique (python?) ! grouper les lignes avec 5% d'erreurs (en fonction de la longueur de la chaîne donc - 50 à 240)

[Totoprosper]

Merci Flodelarab.

En utilisant la commande: sort R1R2_UMI.assembled.Aln.BamTo_[AG]*.fasta | uniq -c > consensus.fasta | grep -v '^>' consensus.fasta > consensus.txt
J'ai effectivement obtenu la listes des séquences avec les occurrences (fréquences) à gauche de chaque séquence.

Maintenant comment se débarrasser dans un premier des fréquences et n'avoir que les séquences seules. Puis dans un deuxième temps, à la place des fréquences, il faudrait séquence_XXXXX avec X=A,C,G, ou T, avec XXXXX indiquant l'origine du fichier dont provient la séquence en question.

Mais les séquences apparemment ne sont pas corrigées. Dans un fichier donné, j'aimerais parcourir et comparer pour chaque groupe de séquences les positions afin de corriger les erreurs de séquençage. Je n'arrive pas faire cela..... une idée?

Merci

[Totoprosper]

Bonjour papajoker,

oui il faut "grouper" dans un fichier donné toutes les séquences ayant au moins 95%. Il faut ensuite corriger les erreurs en se basant sur le fait le séquençage d'une séquence est indépendant de l'autre et il serait impossible voire rare que l'erreur de séquençage arrive à la même position pour deux événement de séquençage différents. Une fois les erreurs corrigées, on prend un représentant de chaque groupe pour mettre mettre dans un nouveau fichier appelé consensus. Donc dans un fichier on pourrait avoir par exemple une centaine de groupes de séquence, voulant signifier qu'on va potentiellement mettre dans le fichier consensus 100 séquences venant de ce fichier... en suite on fait pour les autres fichiers (1024 au total) pour finalement avoir par exemple au final dans le fichier consensus 100*1024 séquences.

Merci

[Flodelarab]

il doit être possible d'avoir un nombre X de caractères différents ("95% uniquement doit être unique")

Entre ATATATATATA et TATATATATA, combien de pourcentage d'erreur ? Et pourquoi ?

[papajoker]

c'est bien possible en python, ici on compare les caractères à la même position:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
a = "ABCDEFGABTDEFGZZUBBAA"
b = "AACXEFGABTDEFGZZUBBAA"
 
nb = round(len(a) *(5/100))  # round() pas bon
print('erreur valide de',nb,'caracteres')
 
diff= [(x,y) for x, y in zip(a, b) if x != y]
print('difference ', diff)
print(len(diff), 'differences')
print(f"est égal ? {len(diff)<=nb}")

le retour donne :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
erreur valide de 1 caracteres
difference  [('B', 'A'), ('D', 'X')]
2 differences
est égal ? False

on arrive donc a avoir un test d'égalité (si j'ai bien compris)

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

estEgal = len([x for x, y in zip(a, b) if x != y]) <= round(len(a) *(5/100))

[BufferBob]

salut,

si on prend deux séquences fictives :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
AGCTGGAGGTTCCCCAGCTCTGCAGCTTCATCCTGAAGACCAGCCA
AAGCTGGAGGTTCCCCAGCTCTGCAGCTTCATCCTGAAGACCAGCCA

on considère qu'elles sont identiques à quel pourcentage ?

Edit: en fait c'est la même question que celle de Flodelarab, sous-entendant qu'il faut une métrique de calcul de similarité, comme une distance de Levenshtein par exemple

[papajoker]

bonjour,

je ne comprends pas vos façons de différencier les chaines ... pour moi c'est clairement les caractères à la même position (donc énormément d'erreurs avec l'exemple de BufferBob)
(Donc le premier critère de différence est bien sûr la longueur ?)

C'est ce que l'on peut voir dans le fichier exemple donné

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
4	CCAAGGATGGTCTTTGATTTCAGGAACAGTGAACCCTTGCTCTACCAATGCTTTCATACGTTTAGCCCAATCTTGGATAGAATATTGCTCTGCAAAGTCAGTAAGGTTCTTAGCAGCAAGAGTCCACCAATCAGAAATGTAGGTGACAAGTGTGCTATAAACC
4	CCAAGGATGGTCTTTGATTTCAGGAACAGTGAACCCTTGCTCTACCAATGCTTTCATACGTTTAGCCCAATCTTGGATAGAATATTGCTCTGCAAAGTCAGTAAGGTTCTTAGCAGCAAGAGTCCACCAATCAGAAATGTAGGTGACAAGTGCGCTATAAACC

ici, uniquement le 11 ème à partir de la fin est différent (T>C)

---

Par contre dans le même fichier, il y a bien un unique exemple très ambiguë ou ici un caractère est "inséré" et nous avons donc 2 chaînes de longueur différentes ! (datas d'entrées fausses ?)

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
1	GGAGGCCGAGGCGGGTGGATCATGAGGTCAGGAGATCAAGACCATCCTGGCTAACACGGTGAAACCCCGTCTCTACTAAAAATATAAAAAATTAGCCAGGCGTCGTGGTGGGCGCCTGTAGTTCCAGCTACTCGTGAGGCTGAGGCAGGAGAATGGTGTGAAC
1	GGAGGCCGAGGCGGGTGGATCATGAGGTCAGGAGATTGAGACCATCCTGGCTAACATGGTAAAACCCCATCTCTACTAAAAATACAAAAAATTTAGCCAGGCATGGTGGCGGGCGCCTGTAGTCCCAGCTACTTGGGAGGCTGAGGCAGGAGAATGGTGTGAAC

[BufferBob]

(...) ici un caractère est "inséré" et nous avons donc 2 chaînes de longueur différentes ! (datas d'entrées fausses ?)

c'est ça, c'est toute la question

il me semble assez intuitivement qu'une mutation sur l'ADN peut très bien se produire par insertion ou retrait d'un codon, et qu'il conviendrait alors de les traiter comme 1 erreur et non pas comme la somme des erreurs dues au désalignement

mais ça n'est que mon intuition, il serait intéressant que le PO tranche ce point

[Flodelarab]

Oui.

Et puis ils viennent d'où ces fichiers ? Moi, j'ai un peu le sentiment qu'on nous demande de re-bricoler le résultat d'une analyse sérieuse de base de données.

Totoprosper peut-il nous expliquer ce qu'il a compris de son travail ?

[disedorgue]

Bonjour,

Personnellement, ce que j'en comprends:

Prenons en exemple, les lignes tagué "1" dans son fichier R1R2_UMI_Merged.assembled.Aln.BamTo_GGAAT.txt.

J'ai retiré le tag et mis chacune des 3 ligne tagué "1" dans un fichier:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
$ cat 1C.txt
GGAGGCCGAGGCGGGTGGATCATGAGGTCAGGAGATCAAGACCATCCTGGCTAACACGGTGAAACCCCGTCTCTACTAAAAATATAAAAAATTAGCCAGGCGTCGTGGTGGGCGCCTGTAGTTCCAGCTACTCGTGAGGCTGAGGCAGGAGAATGGTGTGAAC
$ cat 2C.txt
GGAGGCCGAGGCGGGTGGATCATGAGGTCAGGAGATTGAGACCATCCTGGCTAACATGGTAAAACCCCATCTCTACTAAAAATACAAAAAATTTAGCCAGGCATGGTGGCGGGCGCCTGTAGTCCCAGCTACTTGGGAGGCTGAGGCAGGAGAATGGTGTGAAC
$ cat 3C.txt
GGAGGCCGAGGCGGGTGGATCATGAGGTCAGGAGATCGAGACCATCCTGGCTAACAAGGTGAAACCCCGTCTCTACTAAAAATACAAAAAATTAGCCGGGCGCGGTGGCGGGCGCCTGTAGTCCCAGCTACTCGGGAGGCTGAGGCAGGAGAATGGTGTGAAC

Maintenant, dans son fichier consensus, c'est le 3C.txt qui est considéré, or on peut le voir que c'est le meilleur:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
32
33
34
35
36
$ diff -d -y --suppress-common-lines <(sed 's/./&\n/g' 1C.txt) <(sed 's/./&\n/g' 2C.txt)
C                                                             | T
A                                                             | G
C                                                             | T
G                                                             | A
G                                                             | A
T                                                             | C
                                                              > T
G                                                             | A
C                                                             | G
T                                                             | C
T                                                             | C
C                                                             <
G                                                             <
                                                              > G
                                                              > G
$ diff -d -y --suppress-common-lines <(sed 's/./&\n/g' 1C.txt) <(sed 's/./&\n/g' 3C.txt)
A                                                             | G
C                                                             | A
T                                                             | C
A                                                             | G
T                                                             <
                                                              > G
T                                                             | C
T                                                             | C
T                                                             | G
$ diff -d -y --suppress-common-lines <(sed 's/./&\n/g' 2C.txt) <(sed 's/./&\n/g' 3C.txt)
T                                                             | C
T                                                             | A
A                                                             | G
A                                                             | G
T                                                             <
A                                                             | G
A                                                             | G
T                                                             | C
T                                                             | C

Ici, on voit bien que c'est le 3C.txt qui a le moins de différence avec à la fois le 1C.txt et le 2C.txt

Par contre, avec cette méthode d'analyse, ici on est à environ 5,5% d'erreur pour ce cas (je n'ai pas testé pour les autres)...

Un dernier détail: diff n'est pas commutatif, comparé 2C avec 3C ne donne pas obligatoirement le même nombre de différence que de comparé 3C avec 2C

[disedorgue]

Un petit complément en python:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
>>>from difflib import SequenceMatcher
 
>>>C1="GGAGGCCGAGGCGGGTGGATCATGAGGTCAGGAGATCAAGACCATCCTGGCTAACACGGTGAAACCCCGTCTCTACTAAAAATATAAAAAATTAGCCAGGCGTCGTGGTGGGCGCCTGTAGTTCCAGCTACTCGTGAGGCTGAGGCAGGAGAATGGTGTGAAC"
>>>C2="GGAGGCCGAGGCGGGTGGATCATGAGGTCAGGAGATTGAGACCATCCTGGCTAACATGGTAAAACCCCATCTCTACTAAAAATACAAAAAATTTAGCCAGGCATGGTGGCGGGCGCCTGTAGTCCCAGCTACTTGGGAGGCTGAGGCAGGAGAATGGTGTGAAC"
>>>C3="GGAGGCCGAGGCGGGTGGATCATGAGGTCAGGAGATCGAGACCATCCTGGCTAACAAGGTGAAACCCCGTCTCTACTAAAAATACAAAAAATTAGCCGGGCGCGGTGGCGGGCGCCTGTAGTCCCAGCTACTCGGGAGGCTGAGGCAGGAGAATGGTGTGAAC"
 
>>>s= SequenceMatcher(None,C1,C2)
>>>s.ratio()
0.9235474006116208
 
>>>s= SequenceMatcher(None,C1,C3)
>>>s.ratio()
0.950920245398773
 
>>>s= SequenceMatcher(None,C2,C3)
>>>s.ratio()
0.9480122324159022
 
>>>s= SequenceMatcher(None,C3,C2)
>>>s.ratio()
0.9480122324159022

[Totoprosper]

Merci à toutes/tous pour vos interventions. Les séquences à comparer dans un même fichier doivent être de même taille (longueur). J'ai trouvé et installé le logiciel "gencore" qui a résolu le problème.