Transcript TD 2 Composition en GC de génomes prokaryotes

Version enseignant
Module LV348
TD 2 GC – page 1/6
TD 2 Composition en GC de génomes prokaryotes
Revu dernièrement par Mathilde Carpentier, Cyril Gallut et Joël Pothier
Version du 15 janvier 2014
1
Les données
Nous allons travailler sur le génome d’ Escherichia coli, qui est téléchargeable sur ftp://ftp.
ncbi.nlm.nih.gov/genomes/Bacteria/Escherichia_coli_K_12_substr__MG1655/
NC_000913.fna. Ce fichier est au format fasta. Vous créerez un fichier plus de données plus court
pour tester vos fonctions.
2
2.1
Exercices
Lecture des fichiers
Ecrivez (si cela n’a pas déjà été fait) une fonction python permettant de lire un fichier fasta ne
contenant qu’une séquence et qui retourne la séquence lue.
2.2
Composition globale en GC
Ecrivez une fonction qui prend en entrée un séquence d’ADN et retourne un dictionnaire contenant la composition en A,T,G et C de la séquence. Ce dictionnaire sera donc de la forme : {A :<
f requence des A >, T :< f requence des T >, ...}.
Cette composition vous semble-telle homogène ? Comment le tester statistiquement ? Pourquoi à votre
avis parle-t-on de pourcentage GC (ou AT) ? La composition en GC et AT vous semble-t-elle homogène ?
2.3
Variation du contenu GC le long de la séquence
Le but est d’obtenir le graphique 1(a)
Ecrivez une fonction qui affiche pour une taille de fenêtre, un pas et une séquence donnés la proportion de G+C dans chaque fenêtre. Le prototype de cette fonction est donc
computeSlidingGC(seq,window,step). Elle ne retourne rien mais écrit sur la sortie standard (c’est-à-dire à l’écran) des lignes du type :
0 0.52585
1000 0.52617
2000 0.52645...
Vous pourrez ensuite faire le graphique avec gnuplot (lorsque vous êtes dans gnuplot, tapez plot "<votre
fichier>").
Appliquez cette fonction au génome téléchargé et à une séquence aléatoire de même longueur et
même composition globale en ATGC. Expliquez vos résultats.
c
2013-2014
(by UPMC/Licence de biologie/LV348)
15 janvier 2014
Version enseignant
Module LV348
(a)
TD 2 GC – page 2/6
(b)
F IGURE 1 – Variation du pourcentage GC (a) et du GC skew (b) pour Escherichia coli (NC_000913).
Le pourcentage GC est calulé avec un fenêtre de 100 000 bases et avec un pas de 1000 bases.
2.4 GC Skew
Le but est de maintenant de calculer toujours dans une fenêtre glissante le GC skew qui correspond
G−C
à G+C
. Ecrivez la fonction computeSlidingGCskew qui prend donc en paramêtre la séquence, la
taille de la fenêtre et le pas, et affiche pour chaque fenêtre le GC skew. Faites ensuite le graphique avec
gnuplot pour E. coli et pour une séquence aléatoire. Que remarquez vous ? Expliquez vos résultats.
2.5
Génération de séquences aléatoires
Ecrivez une fonction qui retourne une séquence aléatoire de longueur l donnée en paramètre et
dont la composition en A,T,G,C est aussi donnée en paramètre. Les arguments d’entrée sont donc la
longueur de la séquence et un dictionnaire contenant la composition en A,T,G et C de la séquence à
générer.
Solution :
Toutes les fonctions python sont dans le repertoire code.
c
2013-2014
(by UPMC/Licence de biologie/LV348)
15 janvier 2014
Module LV348
1
Version enseignant
TD 2 GC – page 3/6
#!/sw/bin/python
2
3
4
#import random
from random import *
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#Question 1
#Lecture d’un fichier de seuence au format fasta
def readFasta(nomFi):
#Lecture du contenu du fichier
f=open(nomFi,"r")
lines=f.readlines()
f.close()
#Parsage du du contenu du fichier
seq=""
for l in lines:
if l[0] != ’>’:
#Ne pas oublier d’enlever le retour chariot a la fin des lignes
seq=seq+l[:-1]
return(seq)
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
#Question 2
#Calcul de composition globale en ATGC
def compoGlobale(seq):
resDico={}
for i in seq:
if resDico.has_key(i):
resDico[i]=resDico[i]+1
else:
resDico[i]=0
lg=len(seq)
#
for i in resDico.iterkeys():
for i in resDico.keys():
resDico[i]=float(resDico[i])/lg
return(resDico)
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
#Question 3
#Generation d’une sequence aleatoire de longueur lgSeq
#respectant les proportions donnees dans lesProp
def randomSeq(lgSeq, lesProp):
seq = ""
# lesLettres contiennent tous les symboles possibles
lesLettres=lesProp.keys()
for i in range (0, lgSeq):
#On tire un nombre entre 0.0 et 1.0 au hasard
r = random()
iLettre=0
lim=lesProp[lesLettres[iLettre]]
#On chercher a quelle lettre cela correspond
while iLettre!=len(lesLettres)-1 and r>lim:
iLettre=iLettre+1
#print r,iLettre
lim+=lesProp[lesLettres[iLettre]]
seq=seq+lesLettres[iLettre]
return seq
55
56
c
2013-2014
(by UPMC/Licence de biologie/LV348)
15 janvier 2014
Module LV348
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
Version enseignant
TD 2 GC – page 4/6
#Question 4
#fonction affichant pour une taille de fenetre window, un pas step
#et une sequence seq la proportion de G+C dans chaque fenetre
def computeSlidingGC(seq,window,step):
if len(seq)<window:
print "computeSlidingGC:: Erreur: la sequence doit etre plus longue
que la fenetre"
return
nbGC=0
#Calcul de la premiere fenetre
for i in range(0,window):
if seq[i]==’G’ or seq[i]==’C’:
nbGC=nbGC+1
print "0",float(nbGC)/window
#Calcul des autres fenetres
for i in range(step,len(seq)-window+1,step):
#Pour verification
#GC=0
#for j in range (i, i+window) :
#if seq[j]=="C" or seq[j] == "G" :
#GC+=1
#Ce qu’il faudra soustraire
nbGCDebut=0
for j in range(i-step, i):
if seq[j]==’G’ or seq[j]==’C’:
nbGCDebut=nbGCDebut+1
#Ce qu’il faudra Ajouter
nbGCFin=0
for j in range(i-step+window, i+window):
if seq[j]==’G’ or seq[j]==’C’:
nbGCFin=nbGCFin+1
nbGC=nbGC-nbGCDebut+nbGCFin
#if nbGC != GC:
#print "Pb de calcul du GC !!! nbGC ",nbGC," GC ",GC,i,nbGCDebut,
nbGCFin
print i,float(nbGC)/window
91
92
93
94
95
96
97
98
99
100
101
102
103
104
def computeSlidingGCSimple(seq,window,step):
L=len(seq)
slidingGC=[]
for i in range (0, L-window+1, step):
GC=0
for j in range (i, i+window) :
if seq[j]=="C" or seq[j] == "G" :
GC+=1
print i,float(GC)/window
#slidingGC.append(float(GC)/window)
#return slidingGC
105
106
107
108
109
110
#Question 5
#fonction affichant pour une taille de fenetre window, un pas step
#et une sequence seq la proportion de (G-C)/(G+c) dans chaque fenetre
def computeSlidingGCSkew(seq,window,step):
if len(seq)<window:
c
2013-2014
(by UPMC/Licence de biologie/LV348)
15 janvier 2014
Module LV348
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
Version enseignant
TD 2 GC – page 5/6
print "computeSlidingGC:: Erreur: la sequence doit etre plus longue
que la fenetre"
return
nbG=0
nbC=0
#Calcul de la premiere fenetre
for i in range(0,window):
if seq[i]==’G’:
nbG=nbG+1
if seq[i]==’C’:
nbC=nbC+1
print "0",float(nbG-nbC)/(nbG+nbC)
#Calcul des autres fenetres
for i in range(step,len(seq)-window+1,step):
#Pour verification
#GC=0
#for j in range (i, i+window) :
#if seq[j]=="C" or seq[j] == "G" :
#GC+=1
#Ce qu’il faudra soustraire
nbGDebut=0
nbCDebut=0
for j in range(i-step, i):
if seq[j]==’G’:
nbGDebut=nbGDebut+1
if seq[j]==’C’:
nbCDebut=nbCDebut+1
#Ce qu’il faudra Ajouter
nbGFin=0
nbCFin=0
for j in range(i-step+window, i+window):
if seq[j]==’G’:
nbGFin=nbGFin+1
if seq[j]==’C’:
nbCFin=nbCFin+1
nbG=nbG-nbGDebut+nbGFin
nbC=nbC-nbCDebut+nbCFin
#if nbGC != GC:
#print "Pb de calcul du GC !!! nbGC ",nbGC," GC ",GC,i,nbGCDebut,
nbGCFin
print i,float(nbG-nbC)/(nbG+nbC)
150
151
152
153
154
155
156
157
158
159
#Question 1
#maSeq=readFasta("NC_000913.fna")
#maSeq=readFasta("NC_000964_BSubtilis.fna")
maSeq=readFasta("../data/NC_000949_Bburgdorferi.fna")
#print(maSeq)
print "Le fichier *.fna est lu"
160
161
162
163
164
#Question 2
compoSeq=compoGlobale(maSeq)
#print "Composition de la sequence"
#print compoSeq
c
2013-2014
(by UPMC/Licence de biologie/LV348)
15 janvier 2014
Module LV348
Version enseignant
TD 2 GC – page 6/6
165
166
167
168
169
#########################################
#Les resultats en nombre de nt sont
#{’A’: 1142227, ’C’: 1179553, ’T’: 1140969, ’G’: 1176922}
#Test de Ho : est-ce que la repartition ATGC est homogene ?
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
#chisq.test(c(1142227, 1179553, 1140969, 1176922),p=c(0.25,0.25,0.25,0.25))
#
#
Chi-squared test for given probabilities
#
#data: c(1142227, 1179553, 1140969, 1176922)
#X-squared = 1161.035, df = 3, p-value < 2.2e-16
#Le detail est:
#calcules: res$expected
#[1] 1159918 1159918 1159918 1159918
#sqrt((O-C)^2/C)
# res$residuals
#[1] -16.42603 18.23152 -17.59410 15.78861
#Et donc (O-C)^2/C:
#[1] 269.8145 332.3883 309.5522 249.2802
#####################
186
187
188
189
190
191
#Question 3
print "Generation de la sequence aleatoire"
seqAlea=randomSeq(len(maSeq), compoSeq)
print "Composition de la sequence aleatoire"
print compoGlobale(seqAlea)
192
193
194
195
196
197
print len(maSeq)
#Question 4
print "GC"
computeSlidingGC(maSeq, 500, 50)
#computeSlidingGC(seqAlea, 50000, 5000)
198
199
200
201
202
#Question 5
print "GC skew"
computeSlidingGCSkew(maSeq, 500, 50)
#computeSlidingGCSkew(seqAlea, 50000, 5000)
./code/CorrectionGC.py
c
2013-2014
(by UPMC/Licence de biologie/LV348)
15 janvier 2014