|
0
|
1 #!/usr/bin/env python2.7
|
|
|
2 # -*- coding: utf-8 -*-
|
|
|
3
|
|
|
4 '''
|
|
|
5 Created on sep. 2013
|
|
|
6 @author: rachel legendre
|
|
|
7 @copyright: rachel.legendre@igmors.u-psud.fr
|
|
|
8 @license: GPL v3
|
|
|
9 '''
|
|
|
10
|
|
|
11 from __future__ import division
|
|
|
12 import os, sys, optparse, tempfile, subprocess, re, shutil, commands, urllib, time
|
|
|
13 import itertools
|
|
15
|
14 from math import log10
|
|
0
|
15 from decimal import Decimal
|
|
|
16 from Bio import SeqIO
|
|
|
17 from Bio.Seq import Seq
|
|
|
18 from numpy import arange, std, array, linspace, average
|
|
|
19 #from matplotlib import pyplot as pl
|
|
|
20 import matplotlib
|
|
|
21 matplotlib.use('Agg')
|
|
|
22 import matplotlib.pyplot as pl
|
|
|
23 from matplotlib import font_manager
|
|
|
24 from matplotlib import colors
|
|
|
25 import csv
|
|
17
|
26 from scipy import stats, errstate
|
|
0
|
27 from collections import OrderedDict
|
|
10
|
28 import ribo_functions
|
|
|
29 import HTSeq
|
|
0
|
30 # #libraries for debugg
|
|
13
|
31 #import pdb
|
|
15
|
32 import cPickle
|
|
0
|
33
|
|
|
34 def stop_err(msg):
|
|
|
35 sys.stderr.write("%s\n" % msg)
|
|
|
36 sys.stderr.write("Programme aborted at %s\n" % time.asctime(time.localtime(time.time())))
|
|
|
37 sys.exit()
|
|
|
38
|
|
|
39
|
|
|
40 def init_codon_dict():
|
|
|
41
|
|
|
42 Codon_dict = OrderedDict([('AAA', 0), ('AAC', 0), ('AAG', 0), ('AAT', 0), ('ACA', 0), ('ACC', 0), ('ACG', 0), ('ACT', 0), ('AGA', 0), ('AGC', 0), ('AGG', 0), ('AGT', 0), ('ATA', 0), ('ATC', 0), ('ATG', 0), ('ATT', 0), ('CAA', 0), ('CAC', 0), ('CAG', 0), ('CAT', 0), ('CCA', 0), ('CCC', 0), ('CCG', 0), ('CCT', 0), ('CGA', 0), ('CGC', 0), ('CGG', 0), ('CGT', 0), ('CTA', 0), ('CTC', 0), ('CTG', 0), ('CTT', 0), ('GAA', 0), ('GAC', 0), ('GAG', 0), ('GAT', 0), ('GCA', 0), ('GCC', 0), ('GCG', 0), ('GCT', 0), ('GGA', 0), ('GGC', 0), ('GGG', 0), ('GGT', 0), ('GTA', 0), ('GTC', 0), ('GTG', 0), ('GTT', 0), ('TAA', 0), ('TAC', 0), ('TAG', 0), ('TAT', 0), ('TCA', 0), ('TCC', 0), ('TCG', 0), ('TCT', 0), ('TGA', 0), ('TGC', 0), ('TGG', 0), ('TGT', 0), ('TTA', 0), ('TTC', 0), ('TTG', 0), ('TTT', 0)])
|
|
|
43 return Codon_dict
|
|
|
44
|
|
|
45
|
|
|
46
|
|
|
47 def get_codon_usage(bamfile, GFF, site, kmer, a_pos):
|
|
|
48 '''
|
|
|
49 Read GFF dict and get gene codon usage.
|
|
|
50 Return dict of codons usage
|
|
|
51 '''
|
|
|
52 try:
|
|
|
53 codon = init_codon_dict()
|
|
13
|
54 multi_tag = "XS:i:" ## bowtie Tag
|
|
|
55 tag = "IH:i:1" ## RUM tag
|
|
|
56
|
|
10
|
57 for feature in GFF :
|
|
|
58 if feature.type == 'gene' :
|
|
0
|
59 codon_dict = init_codon_dict()
|
|
10
|
60 chrom = feature.iv.chrom
|
|
|
61 start = feature.iv.start
|
|
|
62 stop = feature.iv.end
|
|
17
|
63 if start+50 < stop-50 :
|
|
13
|
64 region = chrom + ':' + str(start+50) + '-' + str(stop-50)
|
|
17
|
65 # #get all reads in this gene
|
|
|
66 reads = subprocess.check_output(["samtools", "view", bamfile, region])
|
|
|
67 head = subprocess.check_output(["samtools", "view", "-H", bamfile])
|
|
|
68 read_tab = reads.split('\n')
|
|
|
69 for read in read_tab:
|
|
|
70 # # search mapper for eliminate multiple alignements
|
|
|
71 if 'bowtie' in head:
|
|
|
72 multi_tag = "XS:i:"
|
|
|
73 elif 'bwa' in head:
|
|
|
74 multi_tag = "XT:A:R"
|
|
|
75 elif 'TopHat' in head:
|
|
|
76 tag = "NH:i:1"
|
|
|
77 else :
|
|
|
78 stop_err("No PG tag find in "+samfile+". Please use bowtie, bwa or Tophat for mapping")
|
|
|
79
|
|
|
80 if len(read) == 0:
|
|
|
81 continue
|
|
|
82 len_read = len(read.split('\t')[9])
|
|
|
83 # if it's read of good length
|
|
|
84 if len_read == kmer and (tag in read or multi_tag not in read):
|
|
|
85 feat = read.split('\t')
|
|
|
86 seq = feat[9]
|
|
|
87 # if it's a reverse read
|
|
|
88 if feat[1] == '16' :
|
|
|
89 if site == "A" :
|
|
|
90 # #get A-site
|
|
|
91 cod = str(Seq(seq[a_pos-5:a_pos-2]).reverse_complement())
|
|
|
92 elif site == "P" :
|
|
|
93 # #get P-site
|
|
|
94 cod = str(Seq(seq[a_pos-2:a_pos+1]).reverse_complement())
|
|
|
95 else :
|
|
|
96 # #get site-E
|
|
|
97 cod = str(Seq(seq[a_pos+1:a_pos+4]).reverse_complement())
|
|
|
98 # # test if it's a true codon not a CNG codon for example
|
|
|
99 if codon_dict.has_key(cod) :
|
|
|
100 codon_dict[cod] += 1
|
|
|
101 # if it's a forward read
|
|
|
102 elif feat[1] == '0' :
|
|
|
103 if site == "A" :
|
|
|
104 # #get A-site
|
|
|
105 cod = seq[a_pos:a_pos+3]
|
|
|
106 elif site == "P" :
|
|
|
107 # #get P-site
|
|
|
108 cod = seq[a_pos-3:a_pos]
|
|
|
109 else :
|
|
|
110 # #get site-E
|
|
|
111 cod = seq[a_pos-6:a_pos-3]
|
|
|
112 if codon_dict.has_key(cod) :
|
|
|
113 codon_dict[cod] += 1
|
|
|
114 del(read)
|
|
0
|
115 # # add in global dict
|
|
|
116 for cod, count in codon_dict.iteritems() :
|
|
|
117 codon[cod] += count
|
|
17
|
118 if sum(codon.values()) == 0 :
|
|
|
119 stop_err('There are no reads aligning on annotated genes in your GFF file')
|
|
|
120 else :
|
|
|
121 return codon
|
|
0
|
122
|
|
|
123 except Exception, e:
|
|
|
124 stop_err('Error during codon usage calcul: ' + str(e))
|
|
|
125
|
|
|
126
|
|
10
|
127
|
|
|
128
|
|
0
|
129 '''
|
|
|
130 http://pyinsci.blogspot.fr/2009/09/violin-plot-with-matplotlib.html
|
|
|
131 '''
|
|
|
132 def violin_plot(ax, data, pos, bp=False):
|
|
|
133 '''
|
|
|
134 create violin plots on an axis
|
|
|
135 '''
|
|
|
136 dist = max(pos) - min(pos)
|
|
|
137 w = min(0.15 * max(dist, 1.0), 0.5)
|
|
|
138 for d, p in zip(data, pos):
|
|
|
139 k = stats.gaussian_kde(d) # calculates the kernel density
|
|
|
140 m = k.dataset.min() # lower bound of violin
|
|
|
141 M = k.dataset.max() # upper bound of violin
|
|
|
142 x = arange(m, M, (M - m) / 100.) # support for violin
|
|
|
143 v = k.evaluate(x) # violin profile (density curve)
|
|
|
144 v = v / v.max() * w # scaling the violin to the available space
|
|
|
145 ax.fill_betweenx(x, p, v + p, facecolor=color1, alpha=0.3)
|
|
|
146 ax.fill_betweenx(x, p, -v + p, facecolor=color2, alpha=0.3)
|
|
|
147 if bp:
|
|
|
148 ax.boxplot(data, notch=1, positions=pos, vert=1)
|
|
|
149
|
|
|
150
|
|
|
151
|
|
|
152 '''
|
|
|
153 http://log.ooz.ie/2013/02/matplotlib-comparative-histogram-recipe.html
|
|
|
154 '''
|
|
|
155 def comphist(x1, x2, orientation='vertical', **kwargs):
|
|
|
156 """Draw a comparative histogram."""
|
|
|
157 # Split keyword args:
|
|
|
158 kwargs1 = {}
|
|
|
159 kwargs2 = {}
|
|
|
160 kwcommon = {}
|
|
|
161 for arg in kwargs:
|
|
|
162 tgt_arg = arg[:-1]
|
|
|
163 if arg.endswith('1'):
|
|
|
164 arg_dict = kwargs1
|
|
|
165 elif arg.endswith('2'):
|
|
|
166 arg_dict = kwargs2
|
|
|
167 else:
|
|
|
168 arg_dict = kwcommon
|
|
|
169 tgt_arg = arg
|
|
|
170 arg_dict[tgt_arg] = kwargs[arg]
|
|
|
171 kwargs1.update(kwcommon)
|
|
|
172 kwargs2.update(kwcommon)
|
|
|
173
|
|
|
174 fig = pl.figure()
|
|
|
175
|
|
|
176 # Have both histograms share one axis.
|
|
|
177 if orientation == 'vertical':
|
|
|
178 ax1 = pl.subplot(211)
|
|
|
179 ax2 = pl.subplot(212, sharex=ax1)
|
|
|
180 # Flip the ax2 histogram horizontally.
|
|
|
181 ax2.set_ylim(ax1.get_ylim()[::-1])
|
|
|
182 pl.setp(ax1.get_xticklabels(), visible=False)
|
|
|
183 legend_loc = (1, 4)
|
|
|
184 else:
|
|
|
185 ax1 = pl.subplot(122)
|
|
|
186 ax2 = pl.subplot(121, sharey=ax1)
|
|
|
187 # Flip the ax2 histogram vertically.
|
|
|
188 ax2.set_xlim(ax2.get_xlim()[::-1])
|
|
|
189 pl.setp(ax1.get_yticklabels(), visible=False)
|
|
|
190 legend_loc = (1, 2)
|
|
|
191
|
|
|
192 ax1.hist(x1, orientation=orientation, **kwargs1)
|
|
|
193 ax2.hist(x2, orientation=orientation, **kwargs2)
|
|
|
194 ax2.set_ylim(ax1.get_ylim()[::-1])
|
|
|
195 ax1.legend(loc=legend_loc[0])
|
|
|
196 ax2.legend(loc=legend_loc[1])
|
|
|
197 # Tighten up the layout.
|
|
|
198 pl.subplots_adjust(wspace=0.0, hspace=0.0)
|
|
|
199 return fig
|
|
|
200
|
|
|
201
|
|
|
202 def compute_FC_plot(cond1_norm, cond2_norm, cod_name, codon_to_test, dirout):
|
|
|
203
|
|
|
204 FC_tab = []
|
|
|
205 for z, y in zip(cond1_norm.itervalues(), cond2_norm.itervalues()):
|
|
|
206 fc = z - y
|
|
|
207 FC_tab.append(fc)
|
|
|
208 # #codon_to_test = ['TGA','TAG','TAA']
|
|
|
209
|
|
|
210 a = []
|
|
|
211 b = []
|
|
|
212 cod = []
|
|
|
213 for codon in cond1_norm.iterkeys():
|
|
|
214 if codon in codon_to_test :
|
|
|
215 fc = cond1_norm[codon] - cond2_norm[codon]
|
|
|
216 b.append(fc)
|
|
|
217 cod.append(codon)
|
|
|
218 else :
|
|
|
219 fc = cond1_norm[codon] - cond2_norm[codon]
|
|
|
220 a.append(fc)
|
|
|
221
|
|
|
222
|
|
|
223 fig = pl.figure(num=1)
|
|
|
224 comphist(array(a), array(b), label1='All codon', label2=cod_name, color2='green', bins=30, rwidth=1)
|
|
|
225 # pl.show()
|
|
|
226 pl.savefig(dirout + '/hist_codon_fc.png', format="png", dpi=340)
|
|
|
227 pl.clf()
|
|
|
228
|
|
|
229
|
|
|
230 # #violin plot
|
|
|
231 pos = range(2)
|
|
|
232 dat = array([array(a), array(b)])
|
|
|
233 fig = pl.figure()
|
|
|
234 pl.title("Distribution of codons FoldChange between two conditions")
|
|
|
235 ax = fig.add_subplot(1, 1, 1)
|
|
|
236 lab = array(['All codons', cod_name])
|
|
|
237 violin_plot(ax, dat, pos, bp=1)
|
|
|
238 for x, z in zip(dat, pos):
|
|
|
239 ax.plot(z, average(x), color='r', marker='*', markeredgecolor='r')
|
|
|
240 xtickNames = pl.setp(ax, xticklabels=lab)
|
|
|
241 pl.savefig(dirout + '/violinplot_codon.png', format="png", dpi=340)
|
|
|
242 pl.clf()
|
|
|
243
|
|
|
244 # (Fval,pval) = stats.ttest_ind(a, b, axis=0, equal_var=True)
|
|
|
245 (Fval, pval) = stats.mannwhitneyu(a, b)
|
|
|
246 return pval
|
|
|
247
|
|
|
248
|
|
|
249 def get_aa_dict(cond1_norm, cond2_norm):
|
|
|
250
|
|
|
251 # ## create amino acid dictionnary:
|
|
|
252 AA = OrderedDict({})
|
|
|
253 AA['Phe'] = [cond1_norm['TTT'] + cond1_norm['TTC'], cond2_norm['TTT'] + cond2_norm['TTC']]
|
|
|
254 AA['Leu'] = [cond1_norm['TTA'] + cond1_norm['TTG'] + cond1_norm['CTT'] + cond1_norm['CTC'] + cond1_norm['CTA'] + cond1_norm['CTG'], cond2_norm['TTA'] + cond2_norm['TTG'] + cond2_norm['CTT'] + cond2_norm['CTC'] + cond2_norm['CTA'] + cond2_norm['CTG']]
|
|
|
255 AA['Ile'] = [cond1_norm['ATT'] + cond1_norm['ATC'] + cond1_norm['ATA'], cond2_norm['ATT'] + cond2_norm['ATC'] + cond2_norm['ATA']]
|
|
|
256 AA['Met'] = [cond1_norm['ATG'], cond2_norm['ATG']]
|
|
|
257 AA['Val'] = [cond1_norm['GTT'] + cond1_norm['GTC'] + cond1_norm['GTA'] + cond1_norm['GTG'] + cond1_norm['AGT'] + cond1_norm['AGC'], cond2_norm['GTT'] + cond2_norm['GTC'] + cond2_norm['GTA'] + cond2_norm['GTG'] + cond2_norm['AGT'] + cond2_norm['AGC']]
|
|
|
258 AA['Ser'] = [cond1_norm['TCT'] + cond1_norm['TCC'] + cond1_norm['TCA'] + cond1_norm['TCG'], cond2_norm['TCT'] + cond2_norm['TCC'] + cond2_norm['TCA'] + cond2_norm['TCG']]
|
|
|
259 AA['Pro'] = [cond1_norm['CCT'] + cond1_norm['CCC'] + cond1_norm['CCA'] + cond1_norm['CCG'], cond2_norm['CCT'] + cond2_norm['CCC'] + cond2_norm['CCA'] + cond2_norm['CCG']]
|
|
|
260 AA['Thr'] = [cond1_norm['ACT'] + cond1_norm['ACC'] + cond1_norm['ACA'] + cond1_norm['ACG'], cond2_norm['ACT'] + cond2_norm['ACC'] + cond2_norm['ACA'] + cond2_norm['ACG']]
|
|
|
261 AA['Ala'] = [cond1_norm['GCT'] + cond1_norm['GCC'] + cond1_norm['GCA'] + cond1_norm['GCG'], cond2_norm['GCT'] + cond2_norm['GCC'] + cond2_norm['GCA'] + cond2_norm['GCG']]
|
|
|
262 AA['Tyr'] = [cond1_norm['TAT'] + cond1_norm['TAC'], cond2_norm['TAT'] + cond2_norm['TAC']]
|
|
|
263 AA['Stop'] = [cond1_norm['TAA'] + cond1_norm['TAG'] + cond1_norm['TGA'], cond2_norm['TAA'] + cond2_norm['TAG'] + cond2_norm['TGA']]
|
|
|
264 AA['His'] = [cond1_norm['CAT'] + cond1_norm['CAC'], cond2_norm['CAT'] + cond2_norm['CAC']]
|
|
|
265 AA['Gln'] = [cond1_norm['CAA'] + cond1_norm['CAG'], cond2_norm['CAA'] + cond2_norm['CAG']]
|
|
|
266 AA['Asn'] = [cond1_norm['AAT'] + cond1_norm['AAC'], cond2_norm['AAT'] + cond2_norm['AAC']]
|
|
|
267 AA['Lys'] = [cond1_norm['AAA'] + cond1_norm['AAG'], cond2_norm['AAA'] + cond2_norm['AAG']]
|
|
|
268 AA['Asp'] = [cond1_norm['GAT'] + cond1_norm['GAC'], cond2_norm['GAT'] + cond2_norm['GAC']]
|
|
|
269 AA['Glu'] = [cond1_norm['GAA'] + cond1_norm['GAG'], cond2_norm['GAA'] + cond2_norm['GAG']]
|
|
|
270 AA['Cys'] = [cond1_norm['TGT'] + cond1_norm['TGC'], cond2_norm['TGT'] + cond2_norm['TGC']]
|
|
|
271 AA['Trp'] = [cond1_norm['TGG'], cond2_norm['TGG']]
|
|
|
272 AA['Arg'] = [cond1_norm['CGT'] + cond1_norm['CGC'] + cond1_norm['CGA'] + cond1_norm['CGG'] + cond1_norm['AGA'] + cond1_norm['AGG'], cond2_norm['CGT'] + cond2_norm['CGC'] + cond2_norm['CGA'] + cond2_norm['CGG'] + cond2_norm['AGA'] + cond2_norm['AGG']]
|
|
|
273 AA['Gly'] = [cond1_norm['GGT'] + cond1_norm['GGC'] + cond1_norm['GGA'] + cond1_norm['GGG'], cond2_norm['GGT'] + cond2_norm['GGC'] + cond2_norm['GGA'] + cond2_norm['GGG']]
|
|
|
274
|
|
|
275
|
|
|
276 return AA
|
|
|
277
|
|
|
278
|
|
|
279
|
|
|
280 def plot_codon_usage(result, dirout, c1, c2, outfile, color1, color2):
|
|
|
281 '''
|
|
|
282 Take list of dict of codon usage and use matplotlib for do graph
|
|
|
283 '''
|
|
|
284
|
|
|
285 # #if there are replicat
|
|
|
286 if len(result) == 4 :
|
|
|
287 # store each dict in variables to make code more readable
|
|
|
288 cond1_1 = result[0].copy()
|
|
|
289 cond1_2 = result[1].copy()
|
|
|
290 cond2_1 = result[2].copy()
|
|
|
291 cond2_2 = result[3].copy()
|
|
|
292 # get codon order in one of list
|
|
|
293 codon_sorted = sorted(cond1_1.iterkeys(), reverse=False)
|
|
15
|
294 try:
|
|
|
295 # get max of each list
|
|
|
296 sum11 = sum(list(cond1_1.itervalues()))
|
|
|
297 sum12 = sum(list(cond1_2.itervalues()))
|
|
|
298 sum21 = sum(list(cond2_1.itervalues()))
|
|
|
299 sum22 = sum(list(cond2_2.itervalues()))
|
|
|
300 # for each codon, get values and sd in each condition
|
|
|
301 cond1_val = {}
|
|
|
302 cond1 = {}
|
|
|
303 cond2_val = {}
|
|
|
304 cond2 = {}
|
|
|
305 std_cond1 = []
|
|
|
306 std_cond2 = []
|
|
|
307 max_val = [] # # max value for graph
|
|
|
308 for i in codon_sorted:
|
|
17
|
309 # # cond1 = mean of replicats cond1 divided by max
|
|
15
|
310 cond1_val[i] = ((cond1_1[i] / sum11 + cond1_2[i] / sum12) / 2)
|
|
|
311 cond1[i] = ((cond1_1[i] + cond1_2[i]) / 2)
|
|
17
|
312 # # standard deviation = absolute value of difference between replicats of cond1
|
|
15
|
313 std_cond1.append(std(array([(cond1_1[i] * 100 / sum11), (cond1_2[i] * 100 / sum12)])))
|
|
17
|
314 # # cond2 = mean of replicats cond1divided by max
|
|
15
|
315 cond2_val[i] = ((cond2_1[i] / sum21 + cond2_2[i] / sum22) / 2)
|
|
|
316 cond2[i] = ((cond2_1[i] + cond2_2[i]) / 2)
|
|
16
|
317 # # standard deviation = absolute value of difference between replicats of cond2
|
|
15
|
318 std_cond2.append(std(array([((cond2_1[i]) * 100 / sum21), ((cond2_2[i]) * 100 / sum22)])))
|
|
|
319 # # max value for each codon
|
|
|
320 max_val.append(max((cond1_1[i] / sum11 + cond1_2[i] / sum12) / 2, (cond2_1[i] / sum21 + cond2_2[i] / sum22) / 2))
|
|
|
321
|
|
|
322 # for graph design
|
|
|
323 cond1_norm = OrderedDict(sorted(cond1_val.items(), key=lambda t: t[0]))
|
|
|
324 cond1_norm.update ((x, y * 100) for x, y in cond1_norm.items())
|
|
|
325 cond2_norm = OrderedDict(sorted(cond2_val.items(), key=lambda t: t[0]))
|
|
|
326 cond2_norm.update ((x, y * 100) for x, y in cond2_norm.items())
|
|
|
327 max_val = [x * 100 for x in max_val]
|
|
|
328 except ZeroDivisionError:
|
|
|
329 stop_err("Not enough reads to compute the codon occupancy")
|
|
0
|
330
|
|
|
331 AA = get_aa_dict(cond1_norm, cond2_norm)
|
|
|
332 max_valaa = []
|
|
|
333 cond1_aa = []
|
|
|
334 cond2_aa = []
|
|
|
335 aa_name = list(AA.iterkeys())
|
|
|
336 for z in AA.itervalues():
|
|
|
337 cond1_aa.append(z[0])
|
|
|
338 cond2_aa.append(z[1])
|
|
|
339 max_valaa.append(max(z))
|
|
|
340 # # plot amino acid profile :
|
|
15
|
341 fig = pl.figure(figsize=(15,10), num=1)
|
|
0
|
342 width = .50
|
|
|
343 ax = fig.add_subplot(111)
|
|
|
344 ax.xaxis.set_ticks([])
|
|
|
345 ind = arange(21)
|
|
|
346 pl.xlim(0, 21)
|
|
|
347 ax.bar(ind, cond1_aa, width, facecolor=color1, label=c1)
|
|
|
348 ax.bar(ind + width, cond2_aa, width, facecolor=color2, label=c2)
|
|
|
349 for x, y, z in zip(ind, max_valaa, aa_name):
|
|
|
350 ax.text(x + width, y + 0.2, '%s' % z, ha='center', va='bottom', fontsize=14)
|
|
|
351 ax.set_ylabel('Ribosome Occupancy (percent of normalized reads)')
|
|
|
352 ax.set_xlabel('Amino Acid')
|
|
|
353 handles, labels = ax.get_legend_handles_labels()
|
|
|
354 ax.legend(handles, labels)
|
|
|
355 pl.savefig(dirout + '/hist_amino_acid.png', format="png", dpi=340)
|
|
|
356 pl.clf()
|
|
|
357
|
|
|
358
|
|
|
359 # # compute theorical count in COND2
|
|
|
360 sum2 = (sum21 + sum22) / 2
|
|
|
361 cond2_count = []
|
|
|
362 for z in cond1_norm.itervalues() :
|
|
|
363 count = int(z * sum2 / 100)
|
|
|
364 cond2_count.append(count)
|
|
|
365
|
|
|
366 expected = array(cond2_count)
|
|
|
367 observed = array(list(cond2.itervalues()))
|
|
|
368
|
|
|
369 # write result
|
|
|
370 with open(outfile, 'w') as out :
|
|
19
|
371 out.write('Codon,Raw_' + c1 + ',Raw_' + c2 + ',Norm_' + c1 + ',Norm_' + c2 + ',FC(Mut/WT)\n')
|
|
0
|
372 for i in codon_sorted:
|
|
16
|
373 ## if global foldchange is equal to zero
|
|
|
374 if cond1_norm[i] == 0 and cond2_norm[i] == 0:
|
|
19
|
375 out.write(i + ',' + str(cond1[i]) + ',' + str(cond2[i]) + ',' + str(cond1_norm[i]) + ',' + str(cond2_norm[i]) + ',1.0\n')
|
|
16
|
376 elif cond1_norm[i] == 0 :
|
|
19
|
377 out.write(i + ',' + str(cond1[i]) + ',' + str(cond2[i]) + ',' + str(cond1_norm[i]) + ',' + str(cond2_norm[i]) + ',0.0\n')
|
|
16
|
378 else:
|
|
19
|
379 out.write(i + ',' + str(cond1[i]) + ',' + str(cond2[i]) + ',' + str(cond1_norm[i]) + ',' + str(cond2_norm[i]) + ',' + str(cond2_norm[i] / cond1_norm[i]) + '\n')
|
|
17
|
380 with errstate(all='ignore'):
|
|
|
381 chi = stats.chisquare(observed, expected)
|
|
0
|
382 out.write('Khi2 test\n')
|
|
|
383 out.write('T : ' + str(chi[0]) + '; p-value : ' + str(chi[1]) + '\n')
|
|
|
384
|
|
|
385
|
|
|
386
|
|
|
387 # plot result
|
|
15
|
388 fig = pl.figure(figsize=(20,10), num=1)
|
|
13
|
389 width = .40
|
|
0
|
390 ind = arange(len(codon_sorted))
|
|
|
391 ax = fig.add_subplot(111)
|
|
|
392 pl.xlim(0, len(codon_sorted) + 1)
|
|
|
393 ax.spines['right'].set_color('none')
|
|
|
394 ax.spines['top'].set_color('none')
|
|
|
395 ax.xaxis.set_ticks([])
|
|
|
396 ax.spines['left'].set_smart_bounds(True)
|
|
|
397 ax.yaxis.set_ticks_position('left')
|
|
|
398 ax.bar(ind, list(cond1_norm.itervalues()), width, facecolor=color1, yerr=std_cond1, error_kw={'elinewidth':1, 'ecolor':'black'}, label=c1)
|
|
|
399 ax.bar(ind + width, list(cond2_norm.itervalues()), width, yerr=std_cond2, facecolor=color2, error_kw={'elinewidth':1, 'ecolor':'black'}, label=c2)
|
|
|
400 for x, y, z in zip(ind, max_val, codon_sorted):
|
|
|
401 ax.text(x + width, y + 0.2, '%s' % z, ha='center', va='bottom', fontsize=8)
|
|
|
402 ax.set_ylabel('Ribosome Occupancy (percent of normalized reads)')
|
|
|
403 ax.set_xlabel('Codons')
|
|
|
404 handles, labels = ax.get_legend_handles_labels()
|
|
|
405 ax.legend(handles, labels)
|
|
|
406 pl.savefig(dirout + '/hist_codons.png', format="png", dpi=340)
|
|
|
407 pl.clf()
|
|
|
408
|
|
|
409
|
|
|
410 elif len(result) == 2 :
|
|
|
411
|
|
|
412 # store each dict in OrderedDict sorted by key to make code more readable
|
|
|
413 cond1 = result[0]
|
|
|
414 cond2 = result[1]
|
|
|
415 cond1_norm = result[0].copy()
|
|
|
416 cond2_norm = result[1].copy()
|
|
|
417 # pdb.set_trace()
|
|
|
418 # get codon order in one of list
|
|
|
419 codon_sorted = sorted(cond1.iterkeys(), reverse=False)
|
|
15
|
420 try:
|
|
|
421 # get sum of each list
|
|
|
422 sum1 = sum(list(cond1.itervalues()))
|
|
|
423 sum2 = sum(list(cond2.itervalues()))
|
|
|
424 # #Normalize values by sum of each libraries
|
|
|
425 cond1_norm.update ((x, (y / sum1) * 100.0) for x, y in cond1_norm.items())
|
|
|
426 cond2_norm.update((x, (y / sum2) * 100.0) for x, y in cond2_norm.items())
|
|
|
427 except ZeroDivisionError:
|
|
17
|
428 stop_err("Not enough reads to compute the codon occupancy. "+str(sum1)+" and "+str(sum2)+" reads are used for each condition, respectively.\n")
|
|
15
|
429
|
|
0
|
430 # # compute theorical count in COND2
|
|
|
431 cond2_count = []
|
|
|
432 for z in cond1_norm.itervalues() :
|
|
|
433 count = int(z * sum2 / 100.0)
|
|
|
434 cond2_count.append(count)
|
|
|
435
|
|
|
436 expected = array(cond2_count)
|
|
|
437 observed = array(list(cond2.itervalues()))
|
|
|
438
|
|
|
439 AA = get_aa_dict(cond1_norm, cond2_norm)
|
|
|
440
|
|
|
441 max_val = []
|
|
|
442 cond1_aa = []
|
|
|
443 cond2_aa = []
|
|
|
444 aa_name = list(AA.iterkeys())
|
|
|
445 for z in AA.itervalues():
|
|
|
446 cond1_aa.append(z[0])
|
|
|
447 cond2_aa.append(z[1])
|
|
|
448 max_val.append(max(z))
|
|
|
449
|
|
|
450 # # plot amino acid profile :
|
|
15
|
451 fig = pl.figure(figsize=(15,10), num=1)
|
|
10
|
452 width = .45
|
|
0
|
453 ax = fig.add_subplot(111)
|
|
|
454 ind = arange(21)
|
|
|
455 pl.xlim(0, 21)
|
|
|
456 #kwargs = {"hatch":'x'}
|
|
|
457 #ax.bar(ind, cond1_aa, width, facecolor=color1, label=c1, **kwargs)
|
|
|
458 #kwargs = {"hatch":'.'}
|
|
|
459 #ax.bar(ind + width, cond2_aa, width, facecolor=color2, label=c2, **kwargs)
|
|
|
460 ax.bar(ind, cond1_aa, width, facecolor=color1, label=c1)
|
|
|
461 ax.bar(ind + width, cond2_aa, width, facecolor=color2, label=c2)
|
|
|
462 #for x, y, z in zip(ind, max_val, aa_name):
|
|
|
463 # ax.text(x + width, y + 0.2, '%s' % z, ha='center', va='bottom', fontsize=14)
|
|
10
|
464 axis_font = {'size':'10'}
|
|
0
|
465 pl.xticks(ind + width, aa_name,**axis_font)
|
|
|
466 ax.spines['right'].set_visible(False)
|
|
|
467 ax.spines['top'].set_visible(False)
|
|
|
468 ax.yaxis.set_ticks_position('left')
|
|
|
469 ax.xaxis.set_ticks_position('bottom')
|
|
|
470 #ax.xaxis.set_ticks([])
|
|
|
471 ax.set_ylabel('Ribosome Occupancy (percent of normalized reads)',**axis_font)
|
|
|
472 ax.set_xlabel('Amino Acids', **axis_font)
|
|
|
473 handles, labels = ax.get_legend_handles_labels()
|
|
10
|
474 font_prop = font_manager.FontProperties(size=8)
|
|
0
|
475 ax.legend(handles, labels, prop=font_prop)
|
|
|
476 pl.savefig(dirout + '/hist_amino_acid.png', format="png", dpi=340)
|
|
|
477 pl.clf()
|
|
|
478
|
|
|
479 # write result
|
|
|
480 with open(outfile, 'w') as out :
|
|
19
|
481 out.write('Codon,Raw_' + c1 + ',Raw_' + c2 + ',Norm_' + c1 + ',Norm_' + c2 + ',FC(Mut/WT)\n')
|
|
0
|
482 for i in codon_sorted:
|
|
16
|
483 if cond1_norm[i] == 0 and cond2_norm[i] == 0:
|
|
19
|
484 out.write(i + ',' + str(cond1[i]) + ',' + str(cond2[i]) + ',' + str(cond1_norm[i]) + ',' + str(cond2_norm[i]) + ',1.0\n')
|
|
16
|
485 elif cond1_norm[i] == 0 :
|
|
19
|
486 out.write(i + ',' + str(cond1[i]) + ',' + str(cond2[i]) + ',' + str(cond1_norm[i]) + ',' + str(cond2_norm[i]) + ',0.0\n')
|
|
16
|
487 else:
|
|
19
|
488 out.write(i + ',' + str(cond1[i]) + ',' + str(cond2[i]) + ',' + str(cond1_norm[i]) + ',' + str(cond2_norm[i]) + ',' + str(cond2_norm[i] / cond1_norm[i]) + '\n')
|
|
0
|
489 out.write('Khi2 test\n')
|
|
17
|
490 with errstate(all='ignore'):
|
|
|
491 chi = stats.chisquare(observed, expected)
|
|
0
|
492 out.write('T : ' + str(chi[0]) + '; p-value : ' + str(chi[1]) + '\n')
|
|
|
493
|
|
|
494 # # get max value for each codon for histogram
|
|
|
495 max_val = [] # # max value for graph
|
|
|
496 for i in cond1:
|
|
|
497 # # max value for each codon
|
|
|
498 max_val.append(max(cond1_norm[i], cond2_norm[i]))
|
|
|
499
|
|
|
500 # plot result
|
|
15
|
501 fig = pl.figure(figsize=(20,10), num=1)
|
|
10
|
502 #fig = pl.figure(num=1)
|
|
|
503 width = .40
|
|
0
|
504 ind = arange(len(codon_sorted))
|
|
|
505 ax = fig.add_subplot(111)
|
|
|
506 pl.xlim(0, len(codon_sorted) + 1)
|
|
|
507 ax.spines['right'].set_color('none')
|
|
|
508 ax.spines['top'].set_color('none')
|
|
|
509 ax.xaxis.set_ticks([])
|
|
|
510 ax.spines['left'].set_smart_bounds(True)
|
|
|
511 ax.yaxis.set_ticks_position('left')
|
|
|
512 ax.bar(ind, list(cond1_norm.itervalues()), width, facecolor=color1, label=c1)
|
|
|
513 ax.bar(ind + width, list(cond2_norm.itervalues()), width, facecolor=color2, label=c2)
|
|
|
514 for x, y, z in zip(ind, max_val, codon_sorted):
|
|
13
|
515 ax.text(x + width, y + 0.2, '%s' % z, ha='center', va='bottom', fontsize=8)
|
|
0
|
516 ax.set_ylabel('Ribosome Occupancy (percent of normalized reads)')
|
|
|
517 ax.set_xlabel('Codons')
|
|
|
518 handles, labels = ax.get_legend_handles_labels()
|
|
|
519 ax.legend(handles, labels)
|
|
|
520 pl.savefig(dirout + '/hist_codons.png', format="png", dpi=340)
|
|
|
521 pl.clf()
|
|
|
522
|
|
|
523
|
|
|
524 else :
|
|
|
525 stop_err('Error running codon usage plotting : ' + str(e))
|
|
|
526
|
|
|
527
|
|
|
528 return (cond1_norm, cond2_norm, chi[1])
|
|
|
529
|
|
|
530 def write_html_file(html, chi_pval, cond1, cond2):
|
|
|
531 try :
|
|
|
532
|
|
|
533
|
|
|
534 html_str = """
|
|
|
535 <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
|
|
|
536 "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
|
|
|
537
|
|
|
538 <html xmlns="http://www.w3.org/1999/xhtml">
|
|
|
539 <head>
|
|
|
540 <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
|
|
|
541 <link href="/static/june_2007_style/blue/base.css" media="screen" rel="Stylesheet" type="text/css" />
|
|
|
542 </head>
|
|
|
543 <body>
|
|
|
544 <h3>Global visualization</h3>
|
|
|
545 <p>
|
|
15
|
546 <h5>Visualization of density footprint in each codon.</h5><br> If user has selected "Yes" for the replicate option the standard deviation between each replicate is plotted as an error bar in histogram.<br>
|
|
0
|
547 <img border="0" src="hist_codons.png" width="1040"/>
|
|
|
548 </p>
|
|
|
549 <p>
|
|
|
550 <h5>Test for homogeneity distribution between each condition</h5><br>
|
|
|
551 H0 : %s and %s are same distribution <br>
|
|
|
552 Khi2 test p-value: %s<br><br>
|
|
|
553 If p-value less than 0.05, we can reject homogeneity distribution so we can hypothesize that distributions are not the same. Otherwise, we accept H0<br>
|
|
|
554
|
|
|
555 </p>
|
|
|
556 <p>
|
|
|
557 <h5>Visualization of density footprint in each codon groupe by amino acid</h5><br>
|
|
|
558 <img border="0" src="hist_amino_acid.png" width="1040"/>
|
|
|
559 </p>
|
|
|
560 </body>
|
|
|
561 </html> """ % (cond1,cond2,chi_pval)
|
|
|
562
|
|
|
563
|
|
|
564 html_file = open(html, "w")
|
|
|
565 html_file.write(html_str)
|
|
|
566 html_file.close()
|
|
|
567
|
|
|
568 except Exception, e :
|
|
|
569 stop_err('Error during html page creation : ' + str(e))
|
|
|
570
|
|
|
571
|
|
|
572
|
|
|
573
|
|
|
574 def check_codons_list (codons) :
|
|
|
575
|
|
|
576 for codon in codons :
|
|
|
577 if codon not in init_codon_dict().iterkeys() :
|
|
|
578 stop_err('Please to enter a valid codon : ' + codon + ' is not find\n')
|
|
|
579
|
|
|
580
|
|
|
581 def check_index_bam (bamfile) :
|
|
|
582 # #testing indexed bam file
|
|
|
583 if os.path.isfile(bamfile + ".bai") :
|
|
|
584 pass
|
|
|
585 else :
|
|
|
586 cmd = "samtools index %s " % (bamfile)
|
|
|
587 proc = subprocess.Popen(args=cmd, shell=True, stderr=subprocess.PIPE)
|
|
|
588 returncode = proc.wait()
|
|
|
589 # if returncode != 0:
|
|
|
590 # raise Exception
|
|
15
|
591
|
|
|
592 def plot_fc (cond1, cond2, site, dirout):
|
|
|
593
|
|
|
594 fc = cond1.copy()
|
|
|
595
|
|
|
596 for key, value in fc.iteritems():
|
|
17
|
597 if cond1[key] == 0:
|
|
|
598 fc[key] = 1
|
|
|
599 else:
|
|
|
600 fc[key] = cond2[key]/cond1[key]
|
|
15
|
601
|
|
|
602 index = arange(len(fc.keys()))
|
|
|
603 label = fc.keys()
|
|
|
604 label = [w.replace('T','U') for w in label]
|
|
|
605 pl.figure(figsize=(15,10), num=1)
|
|
|
606 ax = pl.subplot(1,1,1)
|
|
|
607 pl.xticks([])
|
|
|
608 pl.scatter(index, fc.values(), color='b')
|
|
|
609 pl.axhline(y=1,color='r')
|
|
|
610 pl.xticks(index, label, rotation=90)
|
|
|
611 pl.ylabel('Foldchange of codon occupancy')
|
|
|
612 ax.yaxis.set_ticks_position('left')
|
|
|
613 ax.xaxis.set_ticks_position('bottom')
|
|
16
|
614 pl.ylim(-1,3)
|
|
15
|
615 pl.title(site+" site")
|
|
|
616 pl.savefig(dirout + '/fc_codons.png', format="png", dpi=340)
|
|
|
617
|
|
0
|
618
|
|
|
619 def __main__():
|
|
10
|
620
|
|
0
|
621
|
|
|
622 # Parse command line options
|
|
|
623 parser = optparse.OptionParser()
|
|
|
624 parser.add_option("-g", "--gff", dest="gff", type="string",
|
|
|
625 help="gff file", metavar="FILE")
|
|
|
626
|
|
|
627 parser.add_option("-1", "--file1", dest="file1", type="string",
|
|
|
628 help="Bam Ribo-Seq alignments cond 1, if rep option, separate files by commas ", metavar="FILE")
|
|
|
629
|
|
|
630 parser.add_option("-2", "--file2", dest="file2", type="string",
|
|
|
631 help="Bam Ribo-Seq alignments cond 2, if rep option, separate files by commas", metavar="FILE")
|
|
|
632
|
|
|
633 parser.add_option("-c", "--cond1", dest="c1", type="string",
|
|
|
634 help="Name for first condition", metavar="STR")
|
|
|
635
|
|
|
636 parser.add_option("-C", "--cond2", dest="c2", type="string",
|
|
|
637 help="Name of second condition", metavar="STR")
|
|
|
638
|
|
13
|
639 parser.add_option("-k", "--kmer", dest="kmer", type="int", default = 28 ,
|
|
10
|
640 help="Length of your phasing reads", metavar="INT")
|
|
0
|
641
|
|
|
642 # parser.add_option("-l", "--list", dest="list_cod", type= "string",
|
|
|
643 # help="list of codons to compare to other", metavar="STR")
|
|
|
644
|
|
|
645 parser.add_option("-o", "--out", dest="outfile", type="string",
|
|
|
646 help="write report to FILE", metavar="FILE")
|
|
|
647
|
|
|
648 parser.add_option("-d", "--dirout", dest="dirout", type="string",
|
|
|
649 help="write report to PNG files", metavar="FILE")
|
|
|
650
|
|
13
|
651 parser.add_option("-a", "--asite", dest="asite", type="int", default = 15 ,
|
|
|
652 help="Off-set from the 5'end of the footprint to the A-site (default is 15)", metavar="INT")
|
|
0
|
653
|
|
13
|
654 parser.add_option("-s", "--site", dest="site", type="string", default = "A" ,
|
|
|
655 help="Script can compute in site A, P or E (default is A-site)", metavar="A|P|E")
|
|
0
|
656
|
|
13
|
657 parser.add_option("-r", "--rep", dest="rep", type="string", default = "no" ,
|
|
0
|
658 help="if replicate or not", metavar="yes|no")
|
|
|
659
|
|
13
|
660 parser.add_option("-x", "--hex_col1", dest="color1", type= "string", default = "SkyBlue" ,
|
|
0
|
661 help="Color for first condition", metavar="STR")
|
|
|
662
|
|
13
|
663 parser.add_option("-X", "--hex_col2", dest="color2", type= "string", default = "Plum" ,
|
|
0
|
664 help="Color for second condition", metavar="STR")
|
|
|
665
|
|
|
666 parser.add_option("-q", "--quiet",
|
|
|
667 action="store_false", dest="verbose", default=True,
|
|
|
668 help="don't print status messages to stdout")
|
|
|
669
|
|
|
670 (options, args) = parser.parse_args()
|
|
|
671 print "Begin codon frequency analysis at", time.asctime(time.localtime(time.time()))
|
|
|
672
|
|
|
673 try:
|
|
|
674 authorized_site = ["A", "P", "E"]
|
|
|
675 if options.site not in authorized_site :
|
|
|
676 stop_err(options.site + ' is not a authorized ribosome site')
|
|
|
677
|
|
|
678 ## Check if colors exist
|
|
|
679 if not colors.is_color_like(options.color1) :
|
|
|
680 stop_err( options.color1+' is not a proper color' )
|
|
|
681 if not colors.is_color_like(options.color2) :
|
|
|
682 stop_err( options.color2+' is not a proper color' )
|
|
|
683
|
|
|
684
|
|
|
685 #### NOT USE IN FINAL VERSION
|
|
|
686 # # get codon list
|
|
|
687 # codons = options.list_cod.upper().split(',')
|
|
|
688 # check_codons_list(codons)
|
|
10
|
689 GFF = HTSeq.GFF_Reader(options.gff)
|
|
0
|
690 # # get html file and directory :
|
|
|
691 (html, html_dir) = options.dirout.split(',')
|
|
17
|
692 if not os.path.exists(html_dir):
|
|
|
693 try:
|
|
|
694 os.mkdir(html_dir)
|
|
|
695 except Exception, e :
|
|
|
696 stop_err('Error running make directory : ' + str(e))
|
|
0
|
697 # #RUN analysis
|
|
|
698 # #If there are replicats
|
|
|
699 if options.rep == "yes" :
|
|
|
700 result = []
|
|
|
701 # split name of each file options by ","
|
|
|
702 cond1 = options.file1.split(',')
|
|
|
703 cond2 = options.file2.split(',')
|
|
|
704 # # calcul for each file
|
|
|
705 for fh in itertools.chain(cond1, cond2):
|
|
|
706 check_index_bam (fh)
|
|
|
707 result.append(get_codon_usage(fh, GFF, options.site, options.kmer, options.asite))
|
|
|
708 (cond1, cond2, chi_pval) = plot_codon_usage(result, html_dir, options.c1, options.c2, options.outfile,options.color1, options.color2)
|
|
|
709 # t_pval = compute_FC_plot(cond1,cond2,codons,html_dir)
|
|
15
|
710 plot_fc (cond1, cond2, options.site, html_dir)
|
|
0
|
711
|
|
|
712 # #If there are no replicat
|
|
|
713 elif options.rep == "no" :
|
|
|
714 result = []
|
|
|
715 # #calcul for each cond
|
|
|
716 for fh in (options.file1, options.file2):
|
|
|
717 check_index_bam (fh)
|
|
|
718 result.append(get_codon_usage(fh, GFF, options.site, options.kmer,options.asite))
|
|
|
719 (cond1, cond2, chi_pval) = plot_codon_usage(result, html_dir, options.c1, options.c2, options.outfile,options.color1, options.color2)
|
|
17
|
720
|
|
0
|
721 # t_pval = compute_FC_plot(cond1,cond2,codons,html_dir)
|
|
15
|
722 plot_fc (cond1, cond2, options.site, html_dir)
|
|
0
|
723 else :
|
|
|
724 sys.stderr.write("Please enter yes or no for --rep option. Programme aborted at %s" % time.asctime(time.localtime(time.time())))
|
|
|
725 sys.exit()
|
|
|
726
|
|
|
727 # write_html_file(html,chi_pval,t_pval,codons,options.c1, options.c2)
|
|
|
728 write_html_file(html, chi_pval, options.c1, options.c2)
|
|
|
729
|
|
|
730 print "Finish codon frequency analysis at", time.asctime(time.localtime(time.time()))
|
|
|
731 except Exception, e:
|
|
|
732 stop_err('Error running codon frequency analysis (main program) : ' + str(e))
|
|
|
733
|
|
|
734
|
|
|
735 if __name__=="__main__":
|
|
|
736 __main__()
|
