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