Discussion :

[Isabella83]

Bonjour à tous,
Je suis débutante en python ( je code en perl habituellement), mais je dois utiliser scipy et donc apprendre python !

Alors, comme exemple de ce que je souhaite faire :

j'ai un fichierA avec des coordonnées :

Code :

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

1
2
chromosome start stop
chrA 1114	1121

qui, en fait, "clusterise" les coordonnées du fichier B ci dessous :

Code :

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

1
2
3
4
5
chromosome start stop count
chrA	1114	1115	2
chrA	1117	1118	2
chrA	1119	1120	1
chrA	1120	1121	1

Et en fait, je voudrais pour chaque cluster dans le fichier A, créer une liste x contenant les start de toutes les coordonnées du fichier B contenu dans le cluster (sachant qu'il faut bien faire attention au chromosome), et une liste y contenant les count de chaque position du fichier B, soit dans l'exemple :

Code :

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

1
2
3
 
X=[1114,1117,1119,11120]
y=[2,2,1,1]

Une fois que j'ai ca, je veux faire une interpolation, mais j'ai déjà le code qui fonctionne, il me manque juste à faire mes listes X et Y en parcourant mes fichiers !

Merci pour vos conseils !

Je ne sais pas du tout comment m'y prendre, avec perl, j'aurai fais ça avec un hash, mais en python ......

[Bulveye]

Salut,

tu peux procéder de la manière suivante:
tu crées une liste vide:
X=[]
with open("input.txt",'r') as f_in:
#ouvre ton fichier texte avec les droits de lecture

with open("output.txt",'w') as f_in:

#ouvre ton fichier texte avec les droits d'ecriture

for line in f_in:

list_line=line.split("ton separateur")

X.append(liste_line[1])

[Isabella83]

J'ai finalement fait comme 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
44
45
46
#!/usr/local/bin/python
 
import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
 
cx = open("cDNA_count.txt", "r")
cluster = open("cluster_on_transcript.txt", "r")
 
 
 
for row in cluster: # parcours du fichier cluster
	x=[]
	y=[]
	ligne = row.strip().split('\t') 
	if len(ligne) == 0 :
		break
	else :
		trans=ligne[0]
		start=ligne[1]
		stop=ligne[2]
 
	for rows in cx : # parcours du fichier avec les cx
		lines = rows.strip().split('\t')
		if len(lines) == 0 :
			break
		else :
			chrsm=lines[0]
			startcx=lines[1]
			stopcx=lines[2]
			count=lines[3]
			if trans == chrsm and start <= startcx and startcx <= stop and start <= stopcx and stopcx <= stop:
				x.append(startcx)
				y.append(count)
 
	tck = interpolate.splrep(x,y)
	xnew = np.arange( np.int(x[0]), np.int(x[len(x)-1]), 0.1)
	ynew = interpolate.splev(xnew,tck,der=0)
	yder = interpolate.splev(xnew,tck,der=1)
	tckder = interpolate.splrep(xnew,yder)
	print interpolate.sproot(tckder)
 
 
 
cx.close()
cluster.close()

Cependant, comment faire pour que ma valeur de x et y changent à chaque ligne parcourue dans le fichier cluster ?

[Isabella83]

J'ai finalement trouvé la solution :

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
#!/usr/local/bin/python
 
import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
 
 
cluster = open("cluster_on_transcript.txt", "r")
 
 
while True :
	row = cluster.readline() # parcours du fichier cluster
	print row
	if row == '' :
		break
	else :
		x=[]
		y=[]
		print x
		print y
		ligne = row.strip().split('\t') 
		trans=ligne[0]
		start=ligne[1]
		stop=ligne[2]
 
	cx = open("cDNA_count.txt", "r")
	for rows in cx : # parcours du fichier avec les cx
		if rows == '' :
			break
		else :
			lines = rows.strip().split('\t')
			chrsm=lines[0]
			startcx=lines[1]
			stopcx=lines[2]
			count=lines[3]
			if trans == chrsm and start <= startcx and startcx <= stop and start <= stopcx and stopcx <= stop:
				x.append(startcx)
				y.append(count)
	print x
	print y
 
 
	tck = interpolate.splrep(x,y)
	xnew = np.arange( np.int(x[0]), np.int(x[len(x)-1]), 0.1)
	ynew = interpolate.splev(xnew,tck,der=0)
	yder = interpolate.splev(xnew,tck,der=1)
	tckder = interpolate.splrep(xnew,yder)
	print interpolate.sproot(tckder)

qu'en pensez vous ?

[__dardanos__]

Salut Isabelle83.
Peut-on avoir de vrais fichiers cluster_on_transcript.txt et cDNA_count.txt pour se faire sa propre idée ?
Ou au moins une partie conséquente de ces fichiers ?

[mont29]

Bonsoir Isabella,

En fait, je vois deux points noirs dans ton code*:

* Tout d’abord, tu fais tes comparaisons de positions sur des chaînes de caractères, et pas des nombres*! Ça doit marcher, mais c’est une solution particulièrement fragile si les données deviennent ne serait-ce que légèrement corrompues/incohérentes, et surtout, c’est horriblement plus lent qu’une comparaison numérique*! Il faudrait donc convertir explicitement tes nombres en int, comme ceci*: trans = int(ligne[0]).

* Pour chaque ligne de ton cluster, tu reparcours et reparses tout ton fichier de “dna count”*! Alors là, si celui-ci devient un tantinet conséquent (genre 100000 lignes ou plus, avec par ex. 10000 clusters)… Pas besoin de faire un dessin (même avec le cache des fichiers de l’OS).

En fait, je lirais/parserais une seule fois le fichier DNA, au début, en stockant ses résultats dans un dict (mapping clé->valeur, équivalent python d’un hash perl, si je ne m’abuse), avec par exemple la position de début comme clé. Tant qu’on reste de l’ordre du million de lignes ça reste raisonnable…

Ou alors, si les enregistrements du fichier cluster et ceux du fichier DNA sont dans le même ordre croissant, ouvrir les deux fichiers au début, et pour chaque cluster, itérer dans les lignes du fichier DNA “juste ce qu’il faut”, c-à-d en s’arrêtant dès que le dernier chromosome du cluster est atteint. C’est probablement la solution la plus performante, même si elle est un peu plus délicate à mettre en œuvre (et présuppose une certaine organisation des données de départ).

[Sve@r]

Salut
Autre facteur d'amélioration: tant qu'à mettre des instructions interruptives qui font hurler les puristes (break/continue), autant utiliser pleinement leur potentiel. Puisque ça interrompt une boucle, ça interrompt alors forcément le if qui y est intégré...

Code python :

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

1
2
3
4
5
while True:
    machin=read()
    if machin == "": break
    ... ici suite du code (else inutile puisque le break a quitté la boucle)...
# while

De plus, j'ai pas vérifié mais il me semble que dans le for rows in cx, le rows ne peut pas être vide...

[PauseKawa]

ben...

Code python :

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

1
2
3
4
5
while True:
    machin=read()
    if machin == "": break
    ... ici suite du code (else inutile puisque le break a quitté la boucle)...
# while

Code :

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

1
2
3
while True:
    # itération

Code :

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

for elem in itérable:

Code :

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

if machin == "": break

Code :

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

1
2
if not machin:
    break  # or continue

Bon, ok ----->[]

[Isabella83]

Alors __dardanos__ :

Voici les fichiers que j'utilise :

cluster_on_transcript.txt

Code :

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

1
2
FBtr0308073	1114	1121
FBtr0070470	2107	2114

et cDNA_count.txt

Code :

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

1
2
3
4
5
6
7
8
FBtr0308073	1114	1115	2
FBtr0308073	1117	1118	2
FBtr0308073	1119	1120	1
FBtr0308073	1120	1121	1
FBtr0070470	2107	2108	3
FBtr0070470	2109	2110	5
FBtr0070470	2111	2112	2
FBtr0070470	2113	2114	1

Mais les fichiers d'origine sont beaucoup plus grands !

Comment puis je optimiser ?
Je ne vois pas comment, en python, stocker mes valeurs dans un dico, pour eviter de parcourir le fichier à chaque fois ?

[Isabella83]

En fait, je lirais/parserais une seule fois le fichier DNA, au début, en stockant ses résultats dans un dict (mapping clé->valeur, équivalent python d’un hash perl, si je ne m’abuse)

Comment puis je faire, car il faut tenir compte du chr + position start et stop ?
Existe t'il un moyen de faire un dico de dico ?

Par exemple , si je veux

dico{chromosome}{start}{stop}=count ?

car tous les champs sont importants dans ce fichier ....

[__dardanos__]

Bonjour Isabella.
Je suis d'accord avec les remarques faites par mont29.

Je propose de faire un traitement avec des tableaux numpy.
Les données sont lues en début de programme avec la fonction loadtxt et converties dans le bon type.
L'extraction des données d'un chromosome avec extract. Il faut utiliser la fonction all pour vérifier les 3 conditions.
Et c'est tout. Après on effectue l'interpolation, puis la mise en forme du résultat (recherche de la valeur maximale par exemple).

Voici un exemple :

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
#!/usr/bin/env python
# -*- coding:Utf-8 -*-
import numpy as np
from scipy import interpolate
 
cluster = np.dtype( [ ('name','|S20'), ('start',np.int), ('stop',np.int) ] )
record  = np.dtype( [ ('name','|S20'), ('start',np.int), ('stop',np.int), ('count',np.int) ] )
 
def extraction_cluster( file_cx, file_cluster ):
    enreg = np.loadtxt( file_cx, dtype=record, skiprows=True)
    resultats = []
    for cx in np.loadtxt( file_cluster, dtype=cluster, skiprows=True) :
        # Extraction des enregistrements
        c1  = cx['name']  == enreg['name']
        c2  = cx['start'] <= enreg['start']
        c3  = cx['stop']  >= enreg['stop']
        res = np.extract( np.all( [c1,c2,c3], axis=0), enreg)
        # Si aucun enregistrement, passe au suivant
        if not res.size: continue
        # Interpolation et recherche des points annulant la dérivée
        tck      = interpolate.splrep( res['start'], res['count'] )
        xnew     = np.linspace( res['start'][0], res['start'][-1], 256)
        yder     = interpolate.splev(xnew,tck,der=1)
        tckder   = interpolate.splrep(xnew,yder)
        racines  = interpolate.sproot(tckder)
        # Valeur maximale
        extremum = interpolate.splev(racines,tck)
        idx      = np.argmax(extremum)
        resultats.append((cx['name'], racines[idx], extremum[idx]))
    return resultats
 
if __name__ == '__main__':
    f_compte  = "cDNA_count.txt"
    f_cluster = "cluster_on_transcript.txt"
    print "Affichage des valeurs maxi"
    for cx, x, y in extraction_cluster( f_compte, f_cluster ):
        print "Chromosome '%s' : %f en %f" % (cx, y, x)

En utilisant l'organisation des données dans le fichier cDNA_count.txt, il est possible d'améliorer l'extraction en évitant de calculer les matrices c1, c2 et c3.

[Isabella83]

merci, j'ai essayé de faire tourner ce script mais j'ai l'erreur suivante :

Code :

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

1
2
3
4
5
6
7
Affichage des valeurs maxi
Traceback (most recent call last):
  File "peak_search.py", line 36, in <module>
    for cx, x, y in extraction_cluster( f_compte, f_cluster ):
  File "peak_search.py", line 12, in extraction_cluster
    for cx in np.loadtxt( file_cluster, dtype=cluster, skiprows=True) :
TypeError: iteration over a 0-d array

[__dardanos__]

Essaie avec ça dans cluster_on_transcript.txt :

Code :

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

1
2
3
chromosome start stop
FBtr0308073 1114    1121
FBtr0070470 2107    2114

et ça dans cDNA_count.txt :

Code :

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

1
2
3
4
5
6
7
8
9
chromosome start stop count
FBtr0308073 1114    1115    2
FBtr0308073 1117    1118    2
FBtr0308073 1119    1120    1
FBtr0308073 1120    1121    1
FBtr0070470 2107    2108    3
FBtr0070470 2109    2110    5
FBtr0070470 2111    2112    2
FBtr0070470 2113    2114    1

[Isabella83]

Ah oui super , ça fonctionne, merci !
Y a t'il une fonction qui existe et qui me permette d'identifier si la valeur de la derivée des point autour des pic potentiels ( identifiés par interpolate.sproot() ) passe d'une valeur >0 à une valeur <0 et ne garder que la valeur du pic correspondant à ce point ?

J'ai fait comme ceci :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
peakPotentiel = interpolate.sproot(tckder)
 
 
	# verifier que la valeur de la derivee passe du positif au negatif autour de la valeur du pic
		for pic in peakPotentiel:
			picStart=pic-1
			picStop = pic+1
			up = interpolate.splev(picStop,tck,der=1)
			down = interpolate.splev(picStart,tck,der=1)
 
			if up <0 and down >0 :
                        print pic

Mais peut etre qu'il y a plus simple ?
Car il se peut qu'il y ait 2 pic par cluster, dans ce cas, comment les identifier avec votre methode car je ne veux pas prendre le max ?

[__dardanos__]

Il faut alors sauver tous les extremum et évaluer la dérivée seconde (np.gradient).
Si elle est négative, le point est un maximum, minimum autrement.
Je n'aurai pas de le temps de t'aider plus aujourd'hui.

[Isabella83]

Merci beaucoup pour ton aide !

Je rencontre un autre problème, qui je pense est dut au fait que certain cluster n'ont que 2 coordonnées (j'ai mis un exemple où ça bloque).

cDNA_count.txt :

Code :

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

1
2
3
4
5
6
7
8
9
10
chr	start	stop	count
FBtr0308073	1114	1115	2
FBtr0308073	1117	1118	2
FBtr0308073	1119	1120	1
FBtr0308073	1120	1121	1
FBtr0070470	2107	2108	3
FBtr0070470	2109	2110	5
FBtr0070470	2111	2112	2
FBtr0070470	2113	2114	1
FBtr0070042	1050	1051	3

cluster_on_transcript.txt

Code :

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

1
2
3
4
chr	start	stop
FBtr0308073	1114	1121
FBtr0070470	2107	2114
FBtr0070042     1050    1051

J'ai l'erreur suivante :

Code :

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

1
2
3
4
5
6
7
8
Traceback (most recent call last):
  File "peak_search.py", line 39, in <module>
    for cx, start,stop,x in extraction_cluster( f_compte, f_cluster ):
  File "peak_search.py", line 24, in extraction_cluster
    tck      = interpolate.splrep( res['start'], res['count'] )
  File "/usr/lib/python2.7/dist-packages/scipy/interpolate/fitpack.py", line 420, in splrep
    raise TypeError('m > k must hold')
TypeError: m > k must hold

J'ai voulu insérer une condition, car j'ai remarqué que dans splrep, il fallait au minimum 3 coordonnées pour interpoler, donc mettre la condition :

Code :

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

if res.size >= 4 :

Mais cela ne fonctionne pas, comment puis je faire ?

[__dardanos__]

Salut.

Les 2 évolutions ont été intégrées :

• Tous les pics d'un chromosome sont sauvegardés (les maximum sont des extremum où la dérivée seconde est négative; il est facile d'évaluer la dérivée avec la spline),
• Les chromosomes ayant moins de 4 enregistrements ne pouvant être traités, un message est affiché et le traitement reprend avec le chromosome suivant (tu avais presque trouvé, il faut mettre un "<" et pas ">=" - l'instruction continue permet de passer à l'élément suivant de la boucle).

Voici le programme :

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
#!/usr/bin/env python
# -*- coding:Utf-8 -*-
import numpy as np
from scipy import interpolate
 
cluster = np.dtype( [ ('name','|S20'), ('start',np.int), ('stop',np.int) ] )
record  = np.dtype( [ ('name','|S20'), ('start',np.int), ('stop',np.int), ('count',np.int) ] )
 
def extraction_cluster( file_cx, file_cluster ):
    enreg = np.loadtxt( file_cx, dtype=record, skiprows=True)
    resultats = []
    for cx in np.loadtxt( file_cluster, dtype=cluster, skiprows=True) :
        # Extraction des enregistrements
        c1  = cx['name']  == enreg['name']
        c2  = cx['start'] <= enreg['start']
        c3  = cx['stop']  >= enreg['stop']
        res = np.extract( np.all( [c1,c2,c3], axis=0), enreg)
        # Si le nombre d'enregistrements est insuffisant, passe au suivant
        if res.size < 4:
            print "Interpolation chromosome '%s' impossible (seulement %d enregistrement(s))" % (cx['name'], res.size)
            continue
        # Interpolation et recherche des points annulant la dérivée
        tck     = interpolate.splrep( res['start'], res['count'] )
        xnew    = np.linspace( res['start'][0], res['start'][-1], 256)
        yder    = interpolate.splev(xnew,tck,der=1)
        tckder  = interpolate.splrep(xnew,yder)
        racines = interpolate.sproot(tckder)
        # Recherche de tous les pics (dérivée seconde négative)
        yder2   = interpolate.splev(racines,tck,der=2)
        x_max   = np.extract( yder2 < 0, racines)
        y_max   = interpolate.splev(x_max,tck)
        # Ajout des pics du chromosome à la liste
        for x, y in zip(x_max,y_max):
            resultats.append((cx['name'], x, y))
    return resultats
 
if __name__ == '__main__':
    f_compte, f_cluster  = "cDNA_count.txt", "cluster_on_transcript.txt"
    print "Affichage des valeurs maxi"
    for cx, x, y in extraction_cluster( f_compte, f_cluster ):
        print "Chromosome '%s' : %f en %f" % (cx, y, x)

[Isabella83]

Bonjour,
Merci beaucoup pour ton aide !

J'ai l'erreur suivante que je ne comprends pas :

Code :

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

1
2
3
4
5
6
7
Affichage des valeurs maxi
Traceback (most recent call last):
  File "peak_search.py", line 40, in <module>
    for cx, x, y in extraction_cluster( f_compte, f_cluster ):
  File "peak_search.py", line 33, in extraction_cluster
    for x, y in zip(x_max,y_max):
TypeError: zip argument #2 must support iteration

[__dardanos__]

C'est difficile de t'aider sans avoir les fichiers d'entrées.
S'il y a un problème sur y_max, c'est que l'interpolation ne s'effectue pas correctement.

Le cas d'une courbe à un extremum qui est un minimum lève une erreur en effet.
Colle ceci après le calcul de x_max :

Code :

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

1
2
3
4
        # cas de la courbe a un seul minimum
        if not x_max.size :
            print "Donnees du chromosome '%s' ne possedent pas de maximum" % cx['name']
            continue

Si c'est pas ça, tu peux identifier le chromosome en question, et regarder s'il n'y a pas un problème sur les données.

[Isabella83]

Mais le fichier que j'utilise est encore le fichier que j'ai mis plus haut ... du coup dès que je lance, j'ai cette erreur, alors je pense que cela bloque dès le premier chr ...

cDNA_count.txt

Code :

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

1
2
3
4
5
6
7
8
9
chr	start	stop	count
FBtr0308073	1114	1115	2
FBtr0308073	1117	1118	2
FBtr0308073	1119	1120	1
FBtr0308073	1120	1121	1
FBtr0070470	2107	2108	3
FBtr0070470	2109	2110	5
FBtr0070470	2111	2112	2
FBtr0070470	2113	2114	1

cluster_on_transcript.txt

Code :

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

1
2
3
chr	start	stop
FBtr0308073	1114	1121
FBtr0070470	2107	2114

J'ai découpé le code, et je pense que cela vient du zip, et surement dut au fait que x_max est de "type" array,

Code :

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

1
2
>>> type(x_max)
<type 'numpy.ndarray'>

alors que
y_max est :

Code :

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

1
2
>>> type(y_max)
<type 'numpy.float64'>