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]

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 ?