Discussion :

[asty69]

Bonjour

j ai un enorme probleme alors j ai un programme en C++ qui calcule des probalité de position de certain element d une chaine.
j ai un programme en python qui fait une HMM et indique les positions de ses probalités sur la chaine

et mon probleme c est que dans ces chaines il y a des gap symolisés par des '-' qui fausse la position reel

et je ne vois pas comment faire en python un code qui replace le tout
en sachant que j'ai la chaine de caractere et les positionner erroné.
si quelqu un a une idée je suis preneur car la c est le black out

merci d avance

[kango]

bonjour,

euh... des détails, des exemples ? j'avoue ne pas avoir compris.

[asty69]

Voila le code C++ qui a partir du ne chaine de caractère (ADN) trouve la probalité de position de certaine séquence de gène

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
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
#include <Seq/AlignedSequenceContainer.h>
#include <Seq/Fasta.h>
#include <Seq/SiteContainerTools.h>
#include <Seq/SequenceContainerTools.h>
 
// The Phyllib library
#include <Phyl/models>
 
#include <Phyl/Newick.h>
#include <Phyl/TreeTemplate.h>
 
#include <Phyl/treelikelihoods>
#include <Phyl/RNonHomogeneousMixedTreeLikelihood.h>
#include <Phyl/SubstitutionModelSet.h>
#include <Phyl/SubstitutionModelSetTools.h>
#include <Phyl/FrequenciesSet.h>
 
#include <Phyl/DRASDRTreeLikelihoodData.h>
 
// The STl library:
#include <cstdlib>
#include <iostream>
 
// The NumCalc library
 
#include <NumCalc/distributions>
 
#include <NumCalc/VectorTools.h>
 
using namespace std;
 
using namespace bpp;
 
int main(int argc, char* argv[])
{
 
  if (argc<7){
    cerr << "string nfseq=argv[1]" << endl;
    cerr << "string nfarbre=argv[2]" << endl;
    cerr << "char* nfout=argv[3]" << endl;
    cerr << "double ka=atof(argv[4])" << endl;
    cerr << "double th=atof(argv[5])" << endl;
    cerr << "double alpha1=atof(argv[6])" << endl;
    cerr << "double beta1=atof(argv[7])" << endl;
    cerr << "double alpha2=atof(argv[8])" << endl;
    cerr << "double beta2=atof(argv[9])" << endl;
    exit(0);
  }
 
  string nfseq=argv[1];
  string nfarbre=argv[2];
  char* nfout=argv[3];
  double ka=atof(argv[4]);
  double th=atof(argv[5]);
  double alpha1=atof(argv[6]);
  double beta1=atof(argv[7]);
  double alpha2=atof(argv[8]);
  double beta2=atof(argv[9]);
 
  const NucleicAlphabet * alphabet = new DNA();
  const RNY * alph=new RNY(*alphabet);
 
  const SequenceTools ST;
  const Fasta f;
 
  SubstitutionModelSetTools SMST;
 
  VectorTools vt;
 
  AlignedSequenceContainer ascRNY(alph);
  AlignedSequenceContainer ascDNA(alphabet);
 
  ascDNA.clear();
  f.read(nfseq,ascDNA);
 
  int nbseq=ascDNA.getNumberOfSequences();
  if (nbseq<1){
    cerr << "Pb lecture " << nfseq << endl;
    exit(0);
  }
 
  int i,j,k,l;
 
  ascRNY.clear();
  for (unsigned int j=0;j<nbseq;j++)
    ascRNY.addSequence(*(ST.RNYslice(ascDNA.getSequence(j))));
 
 
  T92* pT92=new T92(alphabet,ka,th);
  YpR_Sym* pYpR = new YpR_Sym(alph, pT92, 10);
 
  map<string, DiscreteDistribution*> madd;
  madd["rCgT"]=new GammaDiscreteDistribution(7,alpha1,beta1);
  MixedSubstitutionModel* pMSM1=new MixedSubstitutionModel(alph,pYpR,madd);
 
  madd["rCgT"]=new GammaDiscreteDistribution(7,alpha2,beta2);
  MixedSubstitutionModel* pMSM2=new MixedSubstitutionModel(alph,pYpR,madd);
 
  // arbre
 
  TreeTemplate<Node> *pa;
  Newick newick(false,true);
 
  pa = newick.read(nfarbre);
 
  ConstantDistribution cd(1);
 
  map<int, double> freqs;
  SequenceContainerTools::getFrequencies(ascRNY, freqs);
  double t = 0;
  vector<double> rootFreq(alph->getSize());
  for(unsigned int i = 0; i < alph->getSize(); i++) t += freqs[i];
  for(unsigned int i = 0; i < alph->getSize(); i++) rootFreq[i] = freqs[i] / t;
 
 
  FullFrequenciesSet FFR(alph, rootFreq);
  SubstitutionModelSet *  pSMS1=SMST.createHomogeneousModelSet(pMSM1,&FFR,pa);
  SubstitutionModelSet *  pSMS2=SMST.createHomogeneousModelSet(pMSM2,&FFR,pa);
 
  // SubstitutionModel* mm;
  // for ( i=0;i<pMSM->getNumberOfModels();i++){
  //   mm=pMSM->getNModel(i);
  //   vector<string> v=mm->getParameters().getParameterNames();
  //   for ( j=0;j<v.size();j++)
  //     cerr << v[j] << " " << mm->getParameterValue(mm->getParameterNameWithoutNamespace(v[j])) << endl;
 
  // }
 
  int size=ascRNY.getSequence(0).size();
 
  RNonHomogeneousMixedTreeLikelihood* pCTL1=new RNonHomogeneousMixedTreeLikelihood(*pa,ascRNY,pSMS1,&cd);
  RNonHomogeneousMixedTreeLikelihood* pCTL2=new RNonHomogeneousMixedTreeLikelihood(*pa,ascRNY,pSMS2,&cd);
 
  pCTL1->initialize();
  pCTL2->initialize();
 
  Vdouble v1, v2;
  v1=pCTL1->getLogLikelihoodForEachSite();
  v2=pCTL2->getLogLikelihoodForEachSite();
 
  ofstream outf;
  outf.open(nfout);
 
  outf << v1.size() << endl;
  outf << "#1\t#2" << endl;
 
  for (unsigned int i=0;i<size;i++)
    outf << v1[i] << "\t" << v2[i] << endl;
 
  outf.close();
 
  delete pCTL1;
  delete pCTL2;
 
  ////////////////
  // reconstruction
 
  delete alphabet;
  delete alph;
  delete pYpR;
}

Voila le code python qui realise une chaine de Markoff et qui d apres les resultat du C++ donne les positions des elements sur la chaine.

Le probleme c est que dans les Chaines il y a des gap '-'

ex ATCT-----TGTG

et le prog va me dire TTGT commence position 4 et fini 12 il faudrai que ca me corrige en me disan que ca commence position 4 et fini en 8

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
44
45
46
47
48
49
50
import lexique
import partition
import matrice
import commands
import Tkinter
import tkFileDialog
import rpy2
import rpy2.robjects as robjects
import os
 
l=lexique.Lexique(str="1:#1 2:#2")
l.g_inter(1,2,0.001)
l.g_inter(1,1,0.999)
l.g_inter(2,1,0.000008)
l.g_inter(2,2,0.999992)
 
root = Tkinter.Tk()
Fichier_a_traiter = tkFileDialog.askopenfilename(filetypes = [("All", "*"),("Selection du fichier a traiter","*.fa")])
 
 
commands.getstatusoutput('/home/m1sp/bin/decoup_CpG '+ Fichier_a_traiter + ' arbre concat_encodeENm001_ilot.out 3 0.4 2.65 1.5 8.8 1.66')
 
m=matrice.Matrice(fic="concat_encodeENm001_logv.out")
 
m2=matrice.Matrice()
m2.fb(m,l)    #matrice des log-probabilites des etats en chaque position
 
p=partition.Partition()
p.read_Matrice(m2)  # partition des log-probabilites maximales en chaque position
 
# for s in p:
#   if s.num()==[1]:
#       print s.deb(), s.fin(), s.val()/len(s)
 
f=open("concat_encodeENm001_logv.part","w")
f.write(str(p))
f.close()
 
 
# 38307 39910 -0.152090080494
# 293160 295112 -0.115315094923
# 315173 318010 -0.0997288662098
# 438275 440998 -0.127240164583
# 602110 604635 -0.0974600457775
# 666073 666482 -0.649121525171
# 687638 690666 -0.0843770852054
# 1036630 1040301 -0.0908340114119
# 1134827 1136571 -0.258962617417
# 1180215 1181891 -0.123348679174
# 1543759 1546672 -0.102127964673

[kango]

si tu cherches juste à supprimer les "gap":

Code :

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

"ATCT-----TGTG".replace('-','')

[asty69]

non ce que je cherche s est crée un programme qui compte les gap et memorise leur position et qui par rapport a la position des sequence rechercher auparavant donne la position des sequences en deduisant les gap avant et au milieu des positions
7689 5000
4577 3446
2346 1234

mais voila avant au mileu de ses position il y a des GAP
donc avec la connaisance des gap et de leur postion le programme corrige les position

7567 4356
4212 3101
etc

je sais pas si je suis clair, mais tu comprends l idée?

[eyquem]

p=partition.Partition()
for s in p:

p est apparemment une instance de Partition() qui contient plusieurs éléments.
Quelle sorte d’éléments ?



p.read_Matrice(m2)
Que fait la méthode read_Matrice(m2) ?



s.num()
Que fait la méthode num() d’un élément de p ?

[asty69]

p=partition.Partition()
for s in p:

p est apparemment une instance de Partition() qui contient plusieurs éléments.
Quelle sorte d’éléments ?



p.read_Matrice(m2)
Que fait la méthode read_Matrice(m2) ?



s.num()
Que fait la méthode num() d’un élément de p ?

p contient les positions des sequences rechercher.

[asty69]

Bonjour

je sais que je n ai pas été tres claire je vais reformuler tout ca en donnant un exemple simplifié.

j ai un fichier qui contient
ACG- - - - TAGCG-CTA- - - -

via different programme j'obtient les positions des sequences qui m interressent (p)

2 - 5
7 - 9
12 - 15
17 - 20

Donc dans ces positions il compte les gap '-'. bien qu il serai plus simple de les enlever des le debuts, je ne peux pas car ils sont necessaire dans le projet.

apres 3 jours de reflections j ai pensé a un truc
position 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
A C G - - - - T A G C G - C T A - - - -

pd pf
2 - 5
7 - 9
12 - 15
16 - 20

n=0 (n=nombre de gap)
pour i de pd à pf
si pd =lettre
pd=pd

sinon n=n+1
pd= pd - n

sinon
si i= pf
si i= lettre
pf= pf-n
sinon n= n+1
pf= pf-n
sinon si i='-'
n=n+1
pf=pf-n

voila mais une fois que ceci est mis en code il va me manquer le decalage du au gap entre les positions

mes questions c'est comment coder ca sous python et comment faire pour inclure dans le decompte les gap entre les positions

merci de votre aide

[eyquem]

J’ai regardé sur Wikipedia ce que sont les chaînes de Markov. Ouh , c’est chaud !




Ce que je comprends du reste:




p=partition.Partition() crée une instance de classe, mais dans un état naissant qui ne comporte aucune données.
C’est l’activation de la méthode read_Matrice() de l’objet p qui le peuple d’éléments.

Ceux-ci sont eux mêmes des instances d’une certaine classe qui doit être définie et appelée au sein de la méthode read_Matrice() car ces éléments ont des méthodes num() , deb() , fin() , val() .
Ces éléments ne sont donc pas des structures de données Python quelconques, ce sont des p-instances.

L’objet p est lui même une instance qui dispose d’une méthode d’itération puisqu’on peut faire for s in p: pour débiter les p-objets qu’il contient les uns après les autres.

Les p-objets peuvent être vus comme de courtes séquences ADN munies de méthodes. Il me semble logique de penser que la création de chaque p-objets résulte de calculs qui sont faits au sein de l’appel de la méthode read_Matrice() en recourant au code C++ que tu as donné (auquel je ne comprends rien) et que ce code repose sur la théorie des chaînes de Markoff. Mais on n’a rien à faire de l’horlogerie qu’il y a dans read_Matrice() , on se contente de savoir que ça produit des p-objets qui sont des séquences munies de méthodes.







Ceci étant posé, j’espère correctement, le problème est que les valeurs s.deb() et s.fin() d’une p-séquence s sont fausses par rapport à un besoin ultérieur.

Par exemple pour ce qui est d’une p-séquence s1 ’TTGT’ :
dans une chaîne ADN ’ATCT-----TGTG’ le programme donne comme résultat pour print s.deb(),’-’, s.fin() :
3-12
(et non pas 4-12 comme tu l’as écrit).

Tu voudrais donc corriger 3-12 en 3-7.
Je suppose qu’il y a ultérieurement à tout ça l’objectif de réaliser des calculs de comparaison de séquences ADN longues dans différents génômes dans lesquels les gaps ’-’ existant dans un génôme parasitent les calculs car ils constituent des différences non significatives relativement à l’évolution biologique. Donc les séquences de deux génômes sont d’abord débarassées de leurs gaps avant de les comparer. Je suppose que c’est ce qui motive la nécessité d’obtenir des bornes pour les p-séquences qui soient corrigées , pour être utilisables avec les grandes séquences expurgées de leurs gaps.





Le problème que je vois est qu’on ne peut pas corriger en intervenant au niveau des calculs eux-mêmes qui sont faits pour produire ce résultat, puisqu’ils ont lieu au sein de la fonction read_Matrice() , dont je ne sais pas comment elle recourt au code C++ (ce que je suppose qu’elle fait) dont je me sens par ailleurs incapable de modifier l’algorithme pour que les résultats sortent 3-7 et non pas 3-12.

Mais je vois un moyen de s’en sortir en refaisant un traitement pour chaque p-séquence de la p-liste p.
Pour cela il faut connaître la séquence précise de la chaîne ADN longue dans laquelle est faite la recherche de toutes les p-séquences et la séquence précise de chacune des p-séquences.
Je pense que chaque p-séquence possède une méthode sequ() qui fournit cette séquence précise: sek1 = s1.sequ()





Posé ainsi, le problème ne me semble pas trop difficile à résoudre et on pourrait même se donner une chaîne ADN et quelques séquences courtes comme base pour écrire le code, mais je préfère attendre une réponse pour :
- savoir si j’ai bien compris et formulé le problème
- que tu nous fournisses une séquence ADN longue et une petite liste de séquences courtes qui sont recherchées dans la longue.

Après ça, je pense qu’il y aura pas mal de monde qui pourra t’aider à écrire le code de rectification des bornes de chaque p-séquence.


PS
Ah ... pas vu le précédent post avant d’envoyer celui-ci.
Va falloir que je vérifie, à l’avenir.

[eyquem]

Peut être ai-je mal compris l’objectif, mais le problème qui me paraissait compliqué au départ m’est apparu à la suite de mon précédent post assez simple.

Je me suis donné un brin d’ADN et des courtes séquences pour écrire un code de correction des positions de ces séquences dans le brin avec tirets ’-’ en positions dans le même brin sans gaps.

Je l’ai rendu générique et avec un affichage qui s’adapte aux données, pour pouvoir changer le brin et les séquences recherchées.

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
import re
 
 
def correction_de_positions(ADN, chaines):
 
    def afficher(li,brin):
        gf = ( brin + '\n' +\
               b[0]*' ' + brin[b[0]:b[1]] + '\n' +\
               b[0]*' ' + str(b[0]) + (b[1]-b[0]-len(str(b[0])))*' ' + str(b[1]) + '\n'
               for s,b in li )
        print('\n'.join(gf), li, '\n'+len(brin)*'=',  sep='\n')
 
 
    p = [ (seq,re.search('-*'.join(seq),ADN).span()) for seq in chaines ]
    afficher(p,ADN)
 
    iADN_dans_ADNbases = [ iADN for iADN,car in enumerate(ADN) if car!='-' ]
    deek = dict( (iADN,j) for j,iADN in enumerate(iADN_dans_ADNbases) )
 
    p_condensed = [ (seq,(deek[a],deek[b])) for (seq,(a,b)) in p ]
    afficher(p_condensed,ADN.replace('-',''))
 
 
 
ADN = '--ATGTC--TGATC---TA---TGTCG-TAATGAT----TTAGGATAAC'
chaines = [ 'GTCTGA' , 'GTCGTAA' , 'AGGATA' ]
correction_de_positions(ADN, chaines)

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
--ATGTC--TGATC---TA---TGTCG-TAATGAT----TTAGGATAAC
    GTC--TGA
    4       12
 
--ATGTC--TGATC---TA---TGTCG-TAATGAT----TTAGGATAAC
                       GTCG-TAA
                       23      31
 
--ATGTC--TGATC---TA---TGTCG-TAATGAT----TTAGGATAAC
                                         AGGATA
                                         41    47
 
[('GTCTGA', (4, 12)), ('GTCGTAA', (23, 31)), ('AGGATA', (41, 47))]
 
=================================================
ATGTCTGATCTATGTCGTAATGATTTAGGATAAC
  GTCTGA
  2     8
 
ATGTCTGATCTATGTCGTAATGATTTAGGATAAC
             GTCGTAA
             13     20
 
ATGTCTGATCTATGTCGTAATGATTTAGGATAAC
                          AGGATA
                          26    32
 
[('GTCTGA', (2, 8)), ('GTCGTAA', (13, 20)), ('AGGATA', (26, 32))]
 
==================================

Le principe est le suivant

Code :

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

1
2
3
4
5
6
7
8
9
10
11
ADN

- - A T G T C - - T G  A  T  C  -
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14  dans ADN



ADN_bases

A T G T C T G A T C
0 1 2 3 4 5 6 7 8 9 dans ADN_bases

Les index des bases ATGC dans ADN se retrouvent dans la succession suivante dans ADN_bases quand on élimine les index des tirets:

Code :

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

1
2
A T G T C T G  A  T  C
2 3 4 5 6 9 10 11 12 13 dans ADN

Cette succession indique que
la position 2 dans ADN se retrouve en 0 dans ADN_bases,
la position 3 dans ADN se retrouve en 1 dans ADN_bases,
la position 4 dans ADN se retrouve en 2 dans ADN_bases,
la position 5 dans ADN se retrouve en 3 dans ADN_bases,
la position 6 dans ADN se retrouve en 4 dans ADN_bases,
la position 9 dans ADN se retrouve en 5 dans ADN_bases,
etc

Il suffit de dresser un dictionnaire
deek[position dans ADN] = position dans ADN_bases
pour corriger les positions trouvées dans ADN

Maintenant, je ne sais pas si c’est bien ce que tu cherches à obtenir

[eyquem]

Mon code ci-dessus est idiot.

Je construis un dictionnaire
position dans ADN : position dans ADN_bases
pour toutes les “position dans ADN“ qui sont celles de lettres AGCT dans ADN (car évidemment pour les tirets on ne peut pas trouver de valeur “’position dans ADN_bases“ puisque ADN_bases n’a plus de tirets)

Mais on n’a pas besoin de TOUTES les correspondances
position dans ADN : position dans ADN_bases

On n’a besoin de ces correspondances que pour les bornes des chaînes auxquelles on s’intéresse.




Un algorithme qui va au plus direct est préférable:

- pour chaque borne, on retire de sa valeur le nombre de tirets qui précèdent cette borne

- avec Python, calculer un nombre de tirets dans un segment est facile quand on sait que

Code :

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

sum(True True False True False)  # vaut 3

- le nombre de tirets avant une borne b est égal au nombre de tirets avant la précédente borne a + le nombre de tirets entre a et b



Une petite explication sur '-*'.join(seq)
C’est évidemment une RE.
Pour seq égale à ’AGGCTTA’ par exemple, '-*'.join(seq) vaut A-*G-*G-*C-*T-*T-*A
Cette RE permet à la regex construite avec elle de rechercher une séquence 'AGGCTTA' entrelardée de tirets.

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
import re
 
def positions_corrigees(ADN,chaines):
 
    def afficher(p,brin):
        gf = ( brin + '\n' +\
               b1*' ' + brin[b1:b2] + '\n' +\
               b1*' ' + str(b1) + (b2-b1-len(str(b1)))*' ' + str(b2) + '\n'
               for b1,b2 in p )
        print('\n'.join(gf), p, '\n'+len(brin)*'=',  sep='\n')
 
 
    p = dict( (re.search('-*'.join(seq),ADN).span(),seq) for seq in chaines)
    bornes = sum((list(seg) for seg in p),[])
    afficher(p,ADN)
 
    deek = {}
    nbtirets,prec = 0,0
    for x in bornes:
        nbtirets += sum('-'==car for car in ADN[prec:x])
        deek[x] = x - nbtirets
        prec = x
 
    p_condensed = dict( ((deek[a],deek[b]),p[(a,b)]) for (a,b) in p )
    afficher(p_condensed,ADN.replace('-',''))
 
 
 
ADN = '--ATGTC--TGATC---TA---TGTCG-TAATGAT----TTAGGATAAC'
chaines = [ 'GTCTGA' , 'GTCGTAA' , 'AGGATA' ]
positions_corrigees(ADN,chaines)

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
--ATGTC--TGATC---TA---TGTCG-TAATGAT----TTAGGATAAC
    GTC--TGA
    4       12
 
--ATGTC--TGATC---TA---TGTCG-TAATGAT----TTAGGATAAC
                       GTCG-TAA
                       23      31
 
--ATGTC--TGATC---TA---TGTCG-TAATGAT----TTAGGATAAC
                                         AGGATA
                                         41    47
 
{(4, 12): 'GTCTGA', (23, 31): 'GTCGTAA', (41, 47): 'AGGATA'}
 
=================================================
ATGTCTGATCTATGTCGTAATGATTTAGGATAAC
  GTCTGA
  2     8
 
ATGTCTGATCTATGTCGTAATGATTTAGGATAAC
             GTCGTAA
             13     20
 
ATGTCTGATCTATGTCGTAATGATTTAGGATAAC
                          AGGATA
                          26    32
 
{(2, 8): 'GTCTGA', (13, 20): 'GTCGTAA', (26, 32): 'AGGATA'}
 
==================================

[eyquem]

Quand on ne peut plus, on peut encore.



Pas besoin de créer un dictionnaire deek intermédiaire.
Il faut constituer le dictionnaire p_condensed directement au moment de l’itération dans le dictionnaire p :

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
import re
 
def positions_corrigees(ADN,chaines):
 
    def afficher(p,brin):
        gf = ( brin + '\n' +\
               b1*' ' + brin[b1:b2] + '\n' +\
               b1*' ' + str(b1) + (b2-b1-len(str(b1)))*' ' + str(b2) + '\n'
               for b1,b2 in p )
        print('\n'.join(gf), p, '\n'+len(brin)*'=',  sep='\n')
 
 
    p = dict( (re.search('-*'.join(seq),ADN).span(),seq) for seq in chaines)
    afficher(p,ADN)
 
    p_condensed = {}
    nbtirets_b,prec = 0,0
    for (a,b) in sorted(p.keys()):
        nbtirets_a  = nbtirets_b + sum('-'==car for car in ADN[prec:a])
        nbtirets_b  = nbtirets_a + sum('-'==car for car in ADN[   a:b])
        p_condensed[ (a-nbtirets_a, b-nbtirets_b) ] = p[(a,b)]
        prec = b
 
    afficher(p_condensed,ADN.replace('-',''))
 
 
 
ADN = '--ATGTC--TGATC---TA---TGTCG-TAATGAT----TTAGGATAAC'
chaines = [ 'GTCTGA' , 'GTCGTAA' , 'AGGATA' ]
positions_corrigees(ADN,chaines)

[asty69]

merci pour votre aide
voila une chaine adn tres longue en piece jointe

et voila les resultat de p

Code :

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

1
2
3
4
5
6
7
8
9
10
11
38307 39910 -0.152090082288
293160 295112 -0.115315093619
315173 318010 -0.0997288646176
438275 440998 -0.127240162134
602110 604635 -0.0974600454273
666073 666482 -0.649121535526
687638 690666 -0.0843770855616
1036630 1040301 -0.0908340132156
1134827 1136571 -0.25896262529
1180215 1181891 -0.123348679296
1543759 1546672 -0.102127965397

je ne vois pas comment ton code peut me fournir un p'(on va dire) avec les positions corriger.

PS tout le reste de ton resonnant est juste a propos du code pas besoin de comprendre le code en c++ et on ne peut pas intervenir dessus

[asty69]

voila en piece jointe le brin d adn que j ai utilisé
il faut coller adn2 a adn1

[eyquem]

Tes deux fichiers ne sont pas clean, il y a des caractères d’informations avant la chaîne ADN , il y a des fins de ligne \n dans la chaîne et il y a, après comme avant la chaîne, des blancs.

Il m’a fallu prélever les bonnes chaînes dans les deux fichiers que tu as fournis et les accoler avant de pouvoir aller plus loin. Je l’ai fait avec le programme suivant.

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
44
45
46
47
48
49
import re
 
adn = {}
 
for fichier in ('adn1.txt','adn2.txt'):
    with open(fichier) as f:
        uh = f.read()
 
    # examen du contenu brut du fichier
    print('uh =',fichier+'.read()','\nlen(uh) =',len(uh))
    som = 0
    for base in ('A','G','C','T','N','-'):
        n = uh.count(base)
        som += n
        print (repr(base),n)
    print('total =',som)
    print()
 
    # nettoyage des caracteres '\n'
    uht = uh.replace('\n','')
    print("uht = uh.replace('\\n','')",'\n','len(uht) =',len(uht))
 
    # examen de la sequence d'ADN proprement dite   
    x,y = re.search('[AGTCN-]+',uht).span()
    adn[fichier] = uht[x:y]
    print('(x,y) =',repr((x,y)))
    print('len(uht[x:y]) =',len(uht[x:y]))
    som = 0
    for base in ('A','G','C','T','N','-'):
        n = uht[x:y].count(base)
        som += n
        print (repr(base),n)
    print('total =',som)
 
    # verification de l'absence de caracteres autres que AGCTN- dans la chaine ADN
    for u in uht[x:y]:
        if u not in ('A','G','C','T','N','-'):
            print (repr(u),"est un caractere exotique dans le chaine ADN")
    print('\n\n')
 
 
 
with open('adn.txt','w') as g:
    g.write(adn['adn1.txt']+adn['adn2.txt'])
 
 
with open('adn.txt') as f:
    uh = f.read()
    print ("Nouveau fichier adn.txt = adn1.txt + adn2.txt  :  len('adn.txt') =",len(uh))

uh = adn1.txt.read()
len(uh) = 563200
'A' 168405
'G' 100923
'C' 104658
'T' 166107
'N' 248
'-' 6659
total = 547000

uht = uh.replace('\n','')
len(uht) = 552256
(x,y) = (512, 547512)
len(uht[x:y]) = 547000
'A' 168405
'G' 100923
'C' 104658
'T' 166107
'N' 248
'-' 6659
total = 547000



uh = adn2.txt.read()
len(uh) = 614400
'A' 176726
'G' 115962
'C' 112006
'T' 183310
'N' 240
'-' 7367
total = 595611

uht = uh.replace('\n','')
len(uht) = 602486
(x,y) = (512, 596123)
len(uht[x:y]) = 595611
'A' 176726
'G' 115962
'C' 112006
'T' 183310
'N' 240
'-' 7367
total = 595611



Nouveau fichier adn.txt = adn1.txt + adn2.txt : len('adn.txt') = 1142611

En partant de la connaissance des positions dans l’ADN avec gaps ’-’, et non plus de celle de chaînes comme dans mes précédents posts, il faut modifier la fonction positions_corrigees() pour lui faire prendre comme argument une liste de positions et non plus une liste de tuples (positions,séquence).

Le code modifié est le suivant:

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
def positions_corrigees(ADN, positions, nbtirets_b = 0, prec = 0):
    for (a,b) in positions:
        nbtirets_a  = nbtirets_b + sum('-'==car for car in ADN[prec:a])
        nbtirets_b  = nbtirets_a + sum('-'==car for car in ADN[   a:b])
        yield (a-nbtirets_a, b-nbtirets_b)
        prec = b
 
 
 
with open('adn.txt') as f:
    ADN = f.read()
 
res = '''38307 39910 -0.152090082288
293160 295112 -0.115315093619
315173 318010 -0.0997288646176
438275 440998 -0.127240162134
602110 604635 -0.0974600454273
666073 666482 -0.649121535526
687638 690666 -0.0843770855616
1036630 1040301 -0.0908340132156
1134827 1136571 -0.25896262529
1180215 1181891 -0.123348679296
1543759 1546672 -0.102127965397'''
 
positions = [ (int(a),int(b)) for a,b,_ in map(str.split,res.splitlines()) ]
 
pc = positions_corrigees(ADN,positions)
 
print ('Positions initiales'.ljust(26) + 'Positions corrigees',
       'dans ADN avec gaps '.ljust(26) + 'dans ADN sans gaps ','', sep='\n')
print ('\n'.join(repr(avant).ljust(26)+repr(apres) for avant,apres in zip(positions,pc)))

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
Positions initiales       Positions corrigees
dans ADN avec gaps        dans ADN sans gaps 
 
(38307, 39910)            (37903, 39467)
(293160, 295112)          (289381, 291332)
(315173, 318010)          (311196, 314011)
(438275, 440998)          (432673, 435368)
(602110, 604635)          (594872, 597362)
(666073, 666482)          (657975, 658384)
(687638, 690666)          (679204, 682211)
(1036630, 1040301)        (1024178, 1027779)
(1134827, 1136571)        (1120928, 1122636)
(1180215, 1181891)        (1166189, 1167865)
(1543759, 1546672)        (1529733, 1532646)

Je ne vois pas où était la difficulté pour arriver à ça.

[asty69]

Merci pour ton aide
les fichier adn on du etre abime par ma faute car j ai du les coupé pour les compresser et les poster
ce qui me posait probleme c est dans ton poste tu recontruisait la chaine adn et vu la taille genome le faite de reconstruire la chaine pour compter aurait mis des siecles.

Je t'explique le projet car je pense que ca t interesse un peu


http://pbil.univ-lyon1.fr/members/gu...laxy_chr21.maf

voila le code adn du Chromosome 21. pour 3 especes

dans certain alignement on a que 2 espece et d autre une. Celle ci nous interesse pas. je m interesse qu au alignements ou il y a 3sequence

j ai fais un programme en C++ qui lie le fichier maf et recupere les alignement ou il y a 3 sequence. et qui cree un fichier.fa du style
> hg18
ATT -----ATT
>Remac
ATTAGT
>pancro
ATT--AGT

que je fais analyse via un autre programme en C++ afin de calcule des probalité d'apration de certaine sequence qui code pour un gene selon l evolution

Le resultat est donné a un programme python qui fais utilise un chaine de markoff pour situé les positions avec les gap

maintenant j ai les position sans les gap qui me pemettent de comparé avec les banque de donnée

donc voila voilou

il me reste encore 2 3 chose a paufiné pour que ce programme soit parfait.
mais je maitrise pas tres bien python d ou mes soucis

merci encore pour ton aide
je l ai fais en C++ car je n ai pas su le faire en pyton)

[eyquem]

Et donc ? y a-t-il besoin de quelque chose de plus, ou bien c’est bon ?

[asty69]

bein c est presque bon

je suis entrain d essayer de remplacer des parti de code pour le rendre plus"generique"

Code :

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

1
2
3
4
with open('adn.txt') as f:
    ADN = f.read()
 
res = '''s.deb(), s.fin(), s.val()/len(s)'''

le adn.txt est ce que je peux le remplacer par un fichier.fa avec des conditions
lui disant de lire uniquement la ligne correspondant a la sequence hg.18

Code :

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

1
2
3
4
5
6
7
8
 t=(open('concat_encodeENm001.fa','r'))
>>>ADNt=t.readline()
>>> t=ADN.readlin
 
>>> if ADN=='>hg.18':
...     ADN=t.readline()
...     while t[0]!='>':
...             ADN=ADN.append(t)

et la ce le bug


concat_encodeENm001.fa

et fait de la maniere suivante:

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
 
>hg18
TCTTTACATAACAAGGTTGGAGATCTGAATATAAAAATTCCCTATAAGAT
ATAGGCAGAGAACATTCAGCAGATCACATATGCTCATATGTGCCTGAAGC
TCTTATCTGCCTACAGGCTTGGGATGCACAATTGGGAAGAAGGGTGTTGG
GGGTGAGGAGGGGGAATGTGCAGGAACTGGGGAAAGAAGTTCAGATAAGA
CAAGTCCAGGGTAAGGAAGCAGCAGGCACCTGGAAAGGGGATCAATCCAG
TAGTTGATAGAAGTTCTGACACAAACTCCAAATAAGTTGAGAACAAACCT
GAGAACAATTGAATGGGTGCTCCACTTCACACCTTGGGGAACAGTAGAAC
AGCAAGAGCAGATAAACAAGTCCCGTCACATGGAGCAGGGTAGTCCTACA
TCTTCTGCCTGTCTCCTGAGGAACAAGAAATGATACAGAATTGCCCATGT
ATGGGAGAAGGCAGATAGCATGAGCTTACCTGGGGAGTCTGGATAAGGA
 
>panTro2
TCTTTACATAACAAGGTTGGAGATCTGAATATAAAAATTCCCTATAAGAT
ATAGGCAGAGAACATTCAGCAGATCACATATGCTCATATGTGCCTGAAGC
TCTTATCTGCCTACAGGCTTGGGATGCACAATTGGGAAGAAGGGTGTTGG
GGGTGAGGAGGGGGAATGTGCAGGAACTGGGGAAAGAAGTTCAGATAAGA
CAAGTCCAGGGTAAGGAAGCAGCAGGCACCTGGAAAGGGGATCAATCCAG
TAGTTGATAGAAGTTCTGACACAAACTCCAAATAAGTTGAGAACAAACCT
GAGAACAATTGAATGGGTGCTCCACTTCACACCTTGGGGAACAGTAGAAC
AGCAAGAGCAGATAAACAAGTCCCGTCACATGGAGCAGGGTAGTCCTACA
TCTTCTGCCTGTCTCCTGAGGAACAAGAAATGATACAGAATTGCCCATGT
ATGGGAGAAGGCAGATAGCATGAGCTTACCTGGGGAGTCTGGATAAGGA
 
>rheMac2
TCTTTACATAACAAGGTTGGAGATCTGAATATAAAAATTCCCTATAAGAT
ATAGGCAGAGAACATTCAGCAGATCACATATGCTCATATGTGCCTGAAGC
TCTTATCTGCCTACAGGCTTGGGATGCACAATTGGGAAGAAGGGTGTTGG
GGGTGAGGAGGGGGAATGTGCAGGAACTGGGGAAAGAAGTTCAGATAAGA
CAAGTCCAGGGTAAGGAAGCAGCAGGCACCTGGAAAGGGGATCAATCCAG
TAGTTGATAGAAGTTCTGACACAAACTCCAAATAAGTTGAGAACAAACCT
GAGAACAATTGAATGGGTGCTCCACTTCACACCTTGGGGAACAGTAGAAC
AGCAAGAGCAGATAAACAAGTCCCGTCACATGGAGCAGGGTAGTCCTACA
TCTTCTGCCTGTCTCCTGAGGAACAAGAAATGATACAGAATTGCCCATGT
ATGGGAGAAGGCAGATAGCATGAGCTTACCTGGGGAGTCTGGATAAGGA

[eyquem]

Si le fichier se présente vraiment sous l’aspect que tu décris, on ne peut pas parler de
« la ligne correspondant a la sequence hg.18 »puisque celle-ci s’étend sur PLUSIEURS lignes.


Pour obtenir une séquence en fonction de son nom sur la ligne la précédant, on peut faire soit

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
voulu = 'hg18'
 
with open('hg18_panTro2_rheMac2') as f:
    lines = f.read().splitlines()
 
for i in range(len(lines)):
    if voulu in lines[i]:
        a = i+1
        voulu = 'trouve'
    if voulu=='trouve' and lines[i].strip()=='':
        sequence = ''.join(lines[a:i])
        break

soit ainsi:

Code :

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

1
2
3
4
5
6
7
import re
 
with open('hg18_panTro2_rheMac2') as f:
    uh = f.read()
 
voulu = 'hg18'
sequence = re.search(voulu+'\s*((?:[AGCTN-]+\r?\n)+)',uh).group(1).replace('\n','')

Le premier code est vraiment simple. Un développeur, même en n’utilisant pas Python, ne devrait pas avoir de problème pour en trouver un semblable. Et en Python donc, ce devrait être encore plus aisé.

Ne pas y arriver et sortir sans broncher ton deuxième code, qui est d’autant plus une horreur que Python est dit facile à apprendre, ça pose question.

Une personne sachant programmer (par exemple en C) peut, en seulement quelques heures, en apprendre assez sur Python pour l’enseigner comme langage d’introduction à la programmation.

Éric Le Bigot
Linux Magazine. Hors-série Janvier/Février 2009.

[asty69]

Bonjour

je viens d apprendre a programmer pour les besoins de ma these je fini le 21 juin
par moment meme des choses j arrive plus a faire
et mon directoire ne m aide pas bcp car ils sont assez vieux mais admiratif
enfin voila

une derniere chose , J' ai pas fais d'erreur la dedans

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
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
import lexique
import partition
import matrice
import commands
import Tkinter
import tkFileDialog
import rpy2
import rpy2.robjects as robjects
import os
import re
 
l=lexique.Lexique(str="1:#1 2:#2")
l.g_inter(1,2,0.001)
l.g_inter(1,1,0.999)
l.g_inter(2,1,0.000008)
l.g_inter(2,2,0.999992)
 
root = Tkinter.Tk()
Fichier_a_traiter = tkFileDialog.askopenfilename(filetypes = [("All", "*"),("Selection du fichier a traiter","*.fa")])
 
 
commands.getstatusoutput('/home/m1sp/bin/decoup_CpG '+ Fichier_a_traiter + ' arbre Fichier_a_traiter.out 3 0.4 2.65 1.5 8.8 1.66')
 
m=matrice.Matrice(fic="concat_encodeENm001_logv.out")
 
m2=matrice.Matrice()
m2.fb(m,l)    #matrice des log-probabilites des etats en chaque position
 
p=partition.Partition()
p.read_Matrice(m2)  # partition des log-probabilites maximales en chaque position
 
# for s in p:
#   if s.num()==[1]:
#       print s.deb(), s.fin(), s.val()/len(s)
 
def positions_corrigees(ADN, positions, nbtirets_b = 0, prec = 0):
    for (a,b) in positions:
        nbtirets_a  = nbtirets_b + sum('-'==car for car in ADN[prec:a])
        nbtirets_b  = nbtirets_a + sum('-'==car for car in ADN[   a:b])
        yield (a-nbtirets_a, b-nbtirets_b)
        prec = b
 
 
with open('hg18_panTro2_rheMac2') as f:
    ADN = f.read()
 
voulu = 'hg18'
sequence = re.search(voulu+'\s*((?:[AGCTN-]+\r?\n)+)',ADN).group(1).replace('\n','') 
 
res = '''s.deb(), s.fin(), s.val()/len(s)'''
 
positions = [ (int(a),int(b)) for a,b,_ in map(str.split,res.splitlines()) ]
 
pc = positions_corrigees(ADN,positions)
 
print ('Positions initiales'.ljust(26) + 'Positions corrigees',
       'dans ADN avec gaps '.ljust(26) + 'dans ADN sans gaps ','', sep='\n')
print ('\n'.join(repr(avant).ljust(26)+repr(apres) for avant,apres in zip(positions,pc)))
f=open("Fichier_a_traiter_logv.part","w")
f.write(str(p))
f.close()
 
 
 
robjects.r("""source('dess.r')""")

Avant de le lancer sur les genomes
encore mille merci
( je t enverrai une copie avec un fort remerciment en citation de fin de these)

commands.getstatusoutput: est une commande qui vient de import os et qui permet de lancer un programme C via python

et je me serai pas interesser a python sauf que c est un des rare langage ou ont peux faire rapidement est assez facilement de chaine de markoff ( ce qui es tres compliqué)
meme les langage de programmation statistique comme R n 'arrive pas a realisé ce type de HMM quelque sooit les different librairie

Encore mille merci