|
6
|
1 #!/usr/bin/env python2.7.3
|
|
|
2 #-*- coding: utf-8 -*-
|
|
|
3
|
|
|
4 '''
|
|
|
5 Created on Dec. 2013
|
|
|
6 @author: rachel legendre
|
|
|
7 @copyright: rachel.legendre@igmors.u-psud.fr
|
|
|
8 @license: GPL v3
|
|
|
9 '''
|
|
|
10
|
|
|
11 import os, sys, time, optparse, shutil, re, urllib, subprocess, tempfile
|
|
|
12 from urllib import unquote
|
|
|
13 from Bio import SeqIO
|
|
|
14 import csv
|
|
|
15 import pysam
|
|
|
16 import HTSeq
|
|
|
17 #from matplotlib import pyplot as pl
|
|
|
18 import matplotlib
|
|
|
19 matplotlib.use('Agg')
|
|
|
20 import matplotlib.pyplot as pl
|
|
|
21 from numpy import arange
|
|
|
22 from matplotlib import ticker as t
|
|
|
23 from PIL import Image
|
|
10
|
24 import ribo_functions
|
|
|
25 ## suppress matplotlib warnings
|
|
|
26 import warnings
|
|
|
27
|
|
|
28
|
|
6
|
29
|
|
|
30 def stop_err( msg ):
|
|
|
31 sys.stderr.write( "%s\n" % msg )
|
|
|
32 sys.stderr.write( "Programme aborted at %s\n" % time.asctime(time.localtime(time.time())))
|
|
|
33 sys.exit()
|
|
|
34
|
|
|
35
|
|
|
36 def compute_rpkm(length,count_gene,count_tot) :
|
|
|
37
|
|
|
38 try :
|
|
|
39 rpkm = "{0:.4f}".format(count_gene*1000000.0/(count_tot*length))
|
|
|
40 except ArithmeticError :
|
|
|
41 stop_err( 'Illegal division by zero')
|
|
|
42 return float(rpkm)
|
|
|
43
|
|
|
44
|
|
|
45 def find_stop(seq) :
|
|
|
46 '''
|
|
|
47 Find stop codon in a sequence and return position of first nucleotide in stop
|
|
|
48 '''
|
|
|
49 stop_codon = ['TAA','TAG','TGA']
|
|
|
50 for nt in range(0,len(seq)-3,3) :
|
|
|
51 codon = seq[nt:nt+3]
|
|
|
52 if codon in stop_codon :
|
|
|
53 return nt
|
|
|
54
|
|
|
55 return -1
|
|
|
56
|
|
|
57
|
|
|
58 def check_met(seq):
|
|
|
59 '''
|
|
|
60 Boolean function for testing presence or absence of methionine in 5 codons following stop codon
|
|
|
61 '''
|
|
|
62 met = 'ATG'
|
|
|
63 for pos in range(0,15,3) :
|
|
|
64 codon = seq[pos:pos+3]
|
|
|
65 if codon in met :
|
|
|
66 return True
|
|
|
67
|
|
|
68 return False
|
|
|
69
|
|
|
70
|
|
|
71 '''
|
|
|
72 feature.iv is a GenomicInterval object :
|
|
|
73 A GenomicInterval object has the following slots, some of which
|
|
|
74 are calculated from the other:
|
|
|
75
|
|
|
76 chrom: The name of a sequence (i.e., chromosome, contig, or
|
|
|
77 the like).
|
|
|
78 start: The start of the interval. Even on the reverse strand,
|
|
|
79 this is always the smaller of the two values 'start' and 'end'.
|
|
|
80 Note that all positions should be given as 0-based value!
|
|
|
81 end: The end of the interval. Following Python convention for
|
|
|
82 ranges, this in one more than the coordinate of the last base
|
|
|
83 that is considered part of the sequence.
|
|
|
84 strand: The strand, as a single character, '+' or '-'. '.' indicates
|
|
10
|
85 that the strand is irrelevant. (Alternatively, pass a Strand object.)
|
|
6
|
86 length: The length of the interval, i.e., end - start
|
|
|
87 start_d: The "directional start" position. This is the position of the
|
|
|
88 first base of the interval, taking the strand into account. Hence,
|
|
|
89 this is the same as 'start' except when strand == '-', in which
|
|
|
90 case it is end-1.
|
|
|
91 end_d: The "directional end": Usually, the same as 'end', but for
|
|
|
92 strand=='-1', it is start-1.
|
|
|
93
|
|
|
94 '''
|
|
10
|
95 def check_overlapping(gff_reader,chrm,start,stop,strand, name):
|
|
6
|
96
|
|
|
97 #### probleme avec les genes completement inclu...
|
|
|
98
|
|
|
99 iv2 = HTSeq.GenomicInterval(chrm,start,stop,strand)
|
|
|
100 for feature in gff_reader:
|
|
|
101 if feature.type == "gene" :
|
|
10
|
102 if feature.iv.overlaps(iv2) and feature.attr.get('gene_name') != name :
|
|
6
|
103 ## if its a reverse gene, we replace start of extension by start of previous gene
|
|
|
104 if strand == '-' :
|
|
|
105 return (feature.iv.end+3,stop)
|
|
|
106 ## else we replace stop of extension by start of following gene
|
|
|
107 else :
|
|
|
108 return (start,feature.iv.start-3)
|
|
|
109 ## if no overlap are find, we return -1
|
|
|
110 return (start,stop)
|
|
|
111
|
|
|
112
|
|
|
113 def pass_length(start,stop) :
|
|
|
114
|
|
|
115 if (stop-start) > 25 :
|
|
|
116 return True
|
|
|
117 else :
|
|
|
118 return False
|
|
|
119
|
|
|
120
|
|
|
121 def check_homo_coverage(gene,GFF,start,stop, aln) :
|
|
|
122
|
|
|
123 chrom = GFF[gene]['chrom']
|
|
|
124 ## compute nunber of nucleotides in CDS with a coverage equal to zero
|
|
|
125 nt_cds_cov = 0
|
|
|
126 nt_cds_num = 0
|
|
|
127 for i in range(1,GFF[gene]['exon_number']+1) :
|
|
|
128 for z in range(GFF[gene]['exon'][i]['start'],GFF[gene]['exon'][i]['stop']):
|
|
|
129 nt_cds_num += 1
|
|
|
130 if aln.count(chrom,z,z+1) == 0 :
|
|
|
131 nt_cds_cov += 1
|
|
|
132
|
|
|
133 ## compute percent of CDS no covering
|
|
|
134 percent = nt_cds_cov*100/nt_cds_num
|
|
|
135
|
|
|
136 ## compute number of nucleotides with no coverage in extension
|
|
|
137 nt_no_cover = 0
|
|
|
138 for pos in range(start,stop,1) :
|
|
|
139 if aln.count(chrom,pos,pos+1) == 0 :
|
|
|
140 nt_no_cover += 1
|
|
|
141 #print gene, nt_cds_cov, percent, nt_no_cover
|
|
|
142 #percent10 = (stop-start)*50/100
|
|
|
143 if (nt_no_cover*100)/(stop-start) > percent :
|
|
|
144 return False
|
|
|
145 else :
|
|
|
146 return True
|
|
|
147
|
|
|
148 def plot_gene ( aln, gene, start_extension, stop_extension, dirout ) :
|
|
|
149
|
|
10
|
150
|
|
|
151 ### ignore matplotlib warnings for galaxy
|
|
|
152 warnings.filterwarnings('ignore')
|
|
6
|
153 try:
|
|
|
154 strand = gene['strand']
|
|
|
155 len_gene = gene['stop']-gene['start']
|
|
|
156 if strand is "-" :
|
|
|
157 len_ext = gene['stop']-start_extension
|
|
|
158 ##Â coverage in all gene + extension
|
|
|
159 start = start_extension-100
|
|
|
160 stop = gene['stop']+100
|
|
|
161 vector1 = [0]*(stop-start)
|
|
|
162 #for pileupcolumn in aln.pileup( gene['chrom'], start, stop):
|
|
|
163 # vector.append(pileupcolumn.n)
|
|
|
164 for read in aln.fetch(gene['chrom'], start, stop):
|
|
|
165 if read.is_reverse :
|
|
|
166 ## for get footprint in P-site (estimate) we take 3 nt in middle of read
|
|
|
167 #pos = (read.pos+13)-start
|
|
|
168 pos = (read.pos-start) + (read.rlen/2)-1
|
|
|
169 for z in range(pos,(pos+3)) :
|
|
|
170 ## le fetch contient des reads qui chevauchent les 30 nt de la fin du gene mais dont le site P
|
|
|
171 ## se trouve avant notre vector, le z devient negatif et la couverture augmente en fin de vecteur (ce qui est faux)
|
|
|
172 if z > 0 :
|
|
|
173 try :
|
|
|
174 vector1[z] += 1
|
|
|
175 except IndexError :
|
|
|
176 pass
|
|
|
177 vector1.reverse()
|
|
|
178 mean_read = float(sum(vector1))/float(len(vector1))
|
|
|
179 cov = [(x/mean_read) for x in vector1]
|
|
|
180 idx_tot = arange(len(cov))
|
|
|
181 ## coverage in extension
|
|
|
182 start_ext = start_extension-40
|
|
|
183 stop_ext = stop_extension+30
|
|
|
184 vector2 = [0]*(stop_ext-start_ext)
|
|
|
185 #for pileupcolumn in aln.pileup( gene['chrom'], start, stop):
|
|
|
186 # vector.append(pileupcolumn.n)
|
|
|
187 for read in aln.fetch(gene['chrom'], start_extension, stop_ext):
|
|
|
188 if read.is_reverse :
|
|
|
189 ## get footprint in P-site
|
|
|
190 #pos = (read.pos+13)-start_ext
|
|
|
191 pos = (read.pos-start_ext) + (read.rlen/2)-1
|
|
|
192 for z in range(pos,(pos+3)) :
|
|
|
193 if z > 0 :
|
|
|
194 try :
|
|
|
195 vector2[z] += 1
|
|
|
196 except IndexError :
|
|
|
197 pass
|
|
|
198 vector2.reverse()
|
|
|
199 #mean_read = float(sum(vector))/float(len(vector))
|
|
|
200 cov_ext = [(x/mean_read) for x in vector2]
|
|
|
201 _max = max(cov_ext[30::])
|
|
|
202 idx_ext = arange(len(cov_ext))
|
|
|
203
|
|
|
204 else :
|
|
|
205 len_ext = stop_extension-gene['start']
|
|
|
206 start = gene['start']-100
|
|
|
207 stop = stop_extension+100
|
|
|
208 vector = [0]*(stop-start)
|
|
|
209 #for pileupcolumn in aln.pileup( gene['chrom'], start, stop):
|
|
|
210 #vector.append(pileupcolumn.n)
|
|
|
211 for read in aln.fetch(gene['chrom'], start, stop):
|
|
|
212 if not read.is_reverse :
|
|
|
213 ## get footprint in P-site
|
|
|
214 #pos = (read.pos+12)-start
|
|
|
215 pos = (read.pos-start) + (read.rlen/2)-1
|
|
|
216 for z in range(pos,(pos+3)) :
|
|
|
217 if z > 0 :
|
|
|
218 try :
|
|
|
219 vector[z] += 1
|
|
|
220 except IndexError :
|
|
|
221 pass
|
|
|
222 mean_read = float(sum(vector))/float(len(vector))
|
|
|
223 cov = [(x/mean_read) for x in vector]
|
|
|
224 idx_tot = arange(len(cov))
|
|
|
225
|
|
|
226 ## coverage in extension
|
|
|
227 start_ext = gene['stop']-30
|
|
|
228 stop_ext = stop_extension+40
|
|
|
229 vector = [0]*(stop-start_ext)
|
|
|
230 for read in aln.fetch(gene['chrom'], start_ext, stop_extension):
|
|
|
231 if not read.is_reverse :
|
|
|
232 ## get footprint in P-site
|
|
|
233 pos = (read.pos-start_ext) + (read.rlen/2)-1
|
|
|
234 for z in range(pos,(pos+3)) :
|
|
|
235 if z > 0 :
|
|
|
236 try :
|
|
|
237 vector[z] += 1
|
|
|
238 except IndexError :
|
|
|
239 pass
|
|
|
240
|
|
|
241 cov_ext = [(x/mean_read) for x in vector]
|
|
|
242 idx_ext = arange(len(cov_ext))
|
|
|
243 _max = max(cov_ext[30::])
|
|
|
244
|
|
|
245 #### PLOT FIGURE ####
|
|
|
246
|
|
|
247 font = {'family' : 'serif','color': 'grey','weight' : 'normal','size' : 16 }
|
|
|
248
|
|
|
249
|
|
|
250 fig = pl.figure(num=1)
|
|
|
251 ## create a big subplot for common y axis on two last subplot
|
|
|
252 ax = fig.add_subplot(2,1,2)
|
|
|
253 ## hide all spines
|
|
|
254 ax.spines['top'].set_visible(False)
|
|
|
255 ax.spines['bottom'].set_visible(False)
|
|
|
256 ax.spines['left'].set_visible(False)
|
|
|
257 ax.spines['right'].set_visible(False)
|
|
|
258 ax.tick_params(labelcolor='w', top='off', bottom='off', left='off', right='off')
|
|
|
259
|
|
|
260 ## plot gene structure
|
|
|
261 ax1 = fig.add_subplot(3,1,1)
|
|
|
262 ax1.set_title(gene['name'])
|
|
|
263 ## hide all spines
|
|
|
264 ax1.spines['right'].set_color('none')
|
|
|
265 ax1.spines['top'].set_color('none')
|
|
|
266 ax1.spines['left'].set_color('none')
|
|
|
267 ax1.spines['bottom'].set_color('none')
|
|
|
268 ax1.set_ylim(0,5)
|
|
|
269 #set xlim as second plot with 100nt switch
|
|
|
270 ax1.set_xlim(-100,len(idx_tot)-100)
|
|
|
271 ## hide ticks
|
|
|
272 pl.yticks([])
|
|
|
273 pl.xticks([])
|
|
|
274 ## if it's a "simple" gene
|
|
|
275 if gene['exon_number'] == 1 :
|
|
|
276 exon = arange(len_gene)
|
|
|
277 x = [3]*len_gene
|
|
|
278 ax1.plot(exon,x,color='#9B9E9C')
|
|
|
279 ax1.fill_between(exon,2,x,color='#9B9E9C')
|
|
|
280 ## if it's a gene with introns
|
|
|
281 else :
|
|
|
282 ## plot a line representing intron
|
|
|
283 intron = arange(len_gene)
|
|
|
284 y = [2.5]*len_gene
|
|
|
285 ax1.plot(intron,y,color='#9B9E9C')
|
|
|
286
|
|
|
287 ## plotting each intron
|
|
|
288 start = gene['start']
|
|
|
289 if strand == '+' :
|
|
|
290 for exon in range(1,gene['exon_number']+1,1) :
|
|
|
291 len_exon = gene['exon'][exon]['stop']-gene['exon'][exon]['start']
|
|
|
292 idx = gene['exon'][exon]['start'] - start
|
|
|
293 exo = arange(idx,idx+len_exon)
|
|
|
294 x = [3]*len_exon
|
|
|
295 ax1.plot(exo,x,color='#9B9E9C')
|
|
|
296 ax1.fill_between(exo,2,x,color='#9B9E9C')
|
|
|
297 else :
|
|
|
298 ## if it's a reverse gene we must reverse exon's position
|
|
|
299 start = gene['start']
|
|
|
300 tabF = [2.5]*len_gene
|
|
|
301 tabR = [2.5]*len_gene
|
|
|
302 for exon in range(1,gene['exon_number']+1,1) :
|
|
|
303 len_exon = gene['exon'][exon]['stop']-gene['exon'][exon]['start']
|
|
|
304 idx = gene['exon'][exon]['start'] - start
|
|
|
305 exo = arange(idx,idx+len_exon)
|
|
|
306 for z in exo :
|
|
|
307 tabF[z] = 3
|
|
|
308 tabR[z] = 2
|
|
|
309 tabF.reverse()
|
|
|
310 tabR.reverse()
|
|
|
311 #pl.ylim(0,5)
|
|
|
312 ax1.plot(tabF,color='#9B9E9C')
|
|
|
313 ax1.plot(tabR,color='#9B9E9C')
|
|
|
314 x = arange(len(tabR))
|
|
|
315 ax1.fill_between(x,tabF,tabR,color='#9B9E9C')
|
|
|
316
|
|
|
317 ## insert arrows (as genome browser representation)
|
|
|
318 yloc = 2.5
|
|
|
319 narrows = 20
|
|
|
320 exonwidth = .8
|
|
|
321 spread = .4 * len_gene / narrows
|
|
|
322 for i in range(narrows):
|
|
|
323 loc = (float(i) * len_gene / narrows)+ (len(idx_tot)/100)*2
|
|
|
324 if strand == '+' :
|
|
|
325 x = [loc - spread, loc, loc - spread]
|
|
|
326 else:
|
|
|
327 x = [loc - spread, loc, loc - spread]
|
|
|
328 y = [yloc - exonwidth / 5, yloc, yloc + exonwidth / 5]
|
|
|
329 ax1.plot(x, y, lw=1.4, color='#676B69')
|
|
|
330
|
|
|
331
|
|
|
332 # plot coverage in all gene + extension
|
|
|
333 ax2 = fig.add_subplot(3,1,2)
|
|
|
334 ## fixe limit to length of coverage for x axis
|
|
|
335 ax2.set_xlim(0,len(idx_tot))
|
|
|
336 ## fixe 4 ticks and associated labels for y axis
|
|
|
337 ax2.set_yticklabels(arange(0,int(max(cov)+20),int((max(cov)+20)/4)))
|
|
|
338 ax2.set_yticks(arange(0,int(max(cov)+20),int((max(cov)+20)/4)))
|
|
|
339 ### add start and stop of gene in axe in place of number
|
|
|
340 ax2.set_xticks([100,(100+len_gene),(100+len_ext)])
|
|
|
341 ax2.set_xticklabels(["{:,}".format(gene['start']),"{:,}".format(gene['stop']),""])
|
|
|
342 ## hide spines and any axis
|
|
|
343 ax2.spines['right'].set_visible(False)
|
|
|
344 ax2.spines['top'].set_visible(False)
|
|
|
345 ax2.yaxis.set_ticks_position('left')
|
|
|
346 ax2.xaxis.set_ticks_position('bottom')
|
|
|
347 ## plot and fill
|
|
|
348 ax2.plot(idx_tot, cov, color="#CC0011")
|
|
|
349 ax2.fill_between(idx_tot, 0,cov,color="#CC0011")
|
|
|
350 ax2.set_title("Genomic coordinates (%s,%s)"% (gene['strand'],gene['chrom']) ,fontdict={'fontsize':'small'})
|
|
|
351
|
|
|
352 # plot zoom coverage in extension
|
|
|
353 ax3 = fig.add_subplot(3,1,3)
|
|
|
354 ## fixe limit to length of coverage for x axis
|
|
|
355 ax3.set_ylim(0,int(_max))
|
|
|
356 # we hide spines
|
|
|
357 ax3.spines['right'].set_visible(False)
|
|
|
358 ax3.spines['top'].set_visible(False)
|
|
|
359 ax3.yaxis.set_ticks_position('left')
|
|
|
360 ax3.xaxis.set_ticks_position('bottom')
|
|
|
361 ##Â add stop and in-frame stop in x axis
|
|
|
362 pl.xticks([30,( stop_extension-start_extension)+30],["stop codon","next in-frame stop"])
|
|
|
363 ## fixe good position for labels
|
|
|
364 (ax3.get_xticklabels()[0]).set_horizontalalignment('center')
|
|
|
365 (ax3.get_xticklabels()[1]).set_horizontalalignment('left')
|
|
|
366 #if max coverage is lower than 2, we have a illegal division by zero
|
|
|
367 if _max > 2 :
|
|
|
368 ax3.set_yticklabels(arange(0,int(_max+1),int((_max+1)/3)))
|
|
|
369 ax3.set_yticks(arange(0,int(_max+1),int((_max+1)/3)))
|
|
|
370 else :
|
|
|
371 ##
|
|
|
372 ax3.set_ylim(0,_max)
|
|
|
373 ax3.ticklabel_format(style='sci', scilimits=(-2,2), axis='y',useOffset=False)
|
|
|
374
|
|
|
375 ax3.plot(idx_ext, cov_ext, color="#CC0011")
|
|
|
376 ax3.fill_between(idx_ext, 0,cov_ext,color="#CC0011")
|
|
|
377
|
|
|
378
|
|
|
379 ## get scale of subplot 3
|
|
|
380 #ax3.text(ax3.get_xlim()[1]-50,ax3.get_ylim()[1]-1, r'50 nt', fontdict=font)
|
|
|
381 #pl.arrow( ax3.get_xlim()[1]-50, ax3.get_ylim()[1]-2, ax3.get_xlim()[1]-50, 0, fc="grey", ec="grey",lw=2)
|
|
|
382
|
|
|
383 ## set common y label
|
|
|
384 ax.set_ylabel('Normalized footprint counts',labelpad=20)
|
|
|
385
|
|
|
386 ## draw and save plot
|
|
|
387 pl.draw()
|
|
|
388 #pl.show()
|
|
|
389 pl.savefig(dirout+".png",format='png', dpi=300)
|
|
|
390 pl.clf()
|
|
|
391
|
|
|
392
|
|
|
393 ## Make thumbnail for html page
|
|
|
394 infile = dirout+'.png'
|
|
|
395 size = 128, 128
|
|
|
396 im = Image.open(infile)
|
|
|
397 im.thumbnail(size, Image.ANTIALIAS)
|
|
|
398 im.save(dirout+"_thumbnail.png", "PNG")
|
|
10
|
399 warnings.resetwarnings()
|
|
6
|
400
|
|
|
401 except Exception, e:
|
|
|
402 stop_err( 'Error during gene plotting : ' + str( e ) )
|
|
|
403
|
|
|
404
|
|
10
|
405 def compute_analysis(bam, GFF, fasta, gff, dirout, extend) :
|
|
6
|
406
|
|
|
407 try:
|
|
|
408 background_file = dirout+"/background_sequence.fasta"
|
|
|
409 file_back = open(background_file, 'w')
|
|
|
410 file_context = open(dirout+"/stop_context.fasta", 'w')
|
|
|
411 file_extension = open(dirout+"/extensions.fasta", 'w')
|
|
|
412 ## Opening Bam file with pysam librarie
|
|
|
413 pysam.index(bam)
|
|
|
414 aln = pysam.Samfile(bam, "rb",header=True, check_header = True)
|
|
|
415 ## Opening fasta file in a dict with BioPython
|
|
|
416 SeqDict = SeqIO.to_dict(SeqIO.parse(open(fasta,"r"),"fasta"))
|
|
|
417
|
|
|
418 ## count total read in bam file
|
|
|
419 cmd = "samtools view -c %s " % (bam)
|
|
|
420 proc = subprocess.Popen( args=cmd, shell=True, stdout = subprocess.PIPE)
|
|
|
421 count_tot = int(proc.stdout.read().rstrip())
|
|
|
422 returncode = proc.wait()
|
|
|
423
|
|
|
424 ## opening a GFF reader for check overlapping
|
|
|
425 gff_reader = HTSeq.GFF_Reader(gff)
|
|
|
426
|
|
|
427 with open(dirout+"/readthrough_result.csv","w") as out :
|
|
|
428 myheader = ['Gene','Name', 'FAIL', 'Stop context','chrom','start extension','stop extension','length extension','RPKM CDS', 'RPKM extension','ratio','Annotation','sequence']
|
|
|
429 writer = csv.writer(out,delimiter='\t')
|
|
|
430 writer.writerow(myheader)
|
|
|
431 for gene in GFF['order'] :
|
|
10
|
432 #print GFF[gene]
|
|
|
433 ## maybe no start position in GTF file so we must to check and replace
|
|
|
434 exon_number = GFF[gene]['exon_number']
|
|
|
435 try : GFF[gene]['start']
|
|
|
436 except :
|
|
|
437 if GFF[gene]['strand'] == '+' :
|
|
|
438 GFF[gene]['start'] = GFF[gene]['exon'][1]['start']
|
|
|
439 else :
|
|
|
440 GFF[gene]['start'] = GFF[gene]['exon'][exon_number]['stop']
|
|
|
441 ## also for stop coordinates
|
|
|
442 try : GFF[gene]['stop']
|
|
|
443 except :
|
|
|
444 if GFF[gene]['strand'] == '+' :
|
|
|
445 GFF[gene]['stop'] = GFF[gene]['exon'][exon_number]['stop']
|
|
|
446
|
|
|
447 else :
|
|
|
448 GFF[gene]['stop'] = GFF[gene]['exon'][1]['start']
|
|
|
449
|
|
|
450
|
|
6
|
451 indic = ""
|
|
|
452 # compute rpkm of CDS :
|
|
|
453 len_cds = GFF[gene]['stop']-GFF[gene]['start']
|
|
|
454 count_cds = 0
|
|
10
|
455 rpkm_cds = 0
|
|
|
456 count_cds = 0
|
|
|
457 try :
|
|
|
458 ### count method of pysam doesn't strand information
|
|
|
459 if GFF[gene]['strand'] == '+' :
|
|
|
460 for track in aln.fetch(GFF[gene]['chrom'],GFF[gene]['start']+12,GFF[gene]['stop']-15) :
|
|
|
461 if not track.is_reverse :
|
|
|
462 count_cds += 1
|
|
|
463 else :
|
|
|
464 for track in aln.fetch(GFF[gene]['chrom'],GFF[gene]['start']+15,GFF[gene]['stop']-12) :
|
|
|
465 if track.is_reverse :
|
|
|
466 count_cds += 1
|
|
|
467 rpkm_cds = compute_rpkm(len_cds,count_cds,count_tot)
|
|
|
468 except ValueError:
|
|
|
469 ## genere warning about gtf coordinates
|
|
|
470 #warnings.warn("Please check coordinates in GFF : maybe stop or start codon are missing" )
|
|
|
471 pass
|
|
6
|
472 ## Only if gene is translate :
|
|
|
473 if rpkm_cds > 0 and count_cds > 128:
|
|
|
474 ## search footprint in UTR3
|
|
|
475 count = 0
|
|
|
476 try :
|
|
|
477 if GFF[gene]['strand'] == '+' :
|
|
|
478 contexte_stop = str(SeqDict[GFF[gene]['chrom']].seq[GFF[gene]['stop']-6:GFF[gene]['stop']+6])
|
|
|
479 #print gene, contexte_stop
|
|
|
480 start_extension = GFF[gene]['stop']
|
|
|
481 stop_extension = GFF[gene]['stop']+90
|
|
|
482 for track in aln.fetch(GFF[gene]['chrom'],start_extension+15,stop_extension) :
|
|
|
483 if not track.is_reverse :
|
|
|
484 count += 1
|
|
|
485
|
|
|
486 ## get sequence of extension
|
|
|
487 seq = str(SeqDict[GFF[gene]['chrom']].seq[start_extension:stop_extension])
|
|
|
488
|
|
|
489 else :
|
|
|
490 contexte_stop = str(SeqDict[GFF[gene]['chrom']].seq[GFF[gene]['start']-7:GFF[gene]['start']+5].reverse_complement())
|
|
|
491 #print gene, contexte_stop
|
|
|
492 start_extension = GFF[gene]['start']-90
|
|
|
493 stop_extension = GFF[gene]['start']
|
|
|
494 for track in aln.fetch(GFF[gene]['chrom'],start_extension,stop_extension-15) :
|
|
|
495 if track.is_reverse :
|
|
|
496 count += 1
|
|
|
497 ## get sequence of extension
|
|
|
498 seq = str(SeqDict[GFF[gene]['chrom']].seq[start_extension:stop_extension-1].reverse_complement())
|
|
10
|
499
|
|
|
500
|
|
6
|
501 ## if we have coverage after cds stop codon
|
|
|
502 if count > 10 :
|
|
|
503 res = find_stop(seq)
|
|
|
504 if res == -1 :
|
|
|
505 '''
|
|
|
506 Write results with no stop but RPF in extension
|
|
|
507 '''
|
|
|
508 ## check if next in-frame codon is far than 90 nt extension :
|
|
|
509 if GFF[gene]['strand'] == '+' :
|
|
10
|
510 pos = check_overlapping(gff_reader,GFF[gene]['chrom'],GFF[gene]['stop']+1,GFF[gene]['stop']+extend,'+',GFF[gene]['name'])
|
|
6
|
511 start_extension = pos[0]-1
|
|
|
512 stop_extension = pos[1]
|
|
10
|
513 #start_extension = GFF[gene]['stop']
|
|
|
514 #stop_extension = GFF[gene]['stop']+extend
|
|
6
|
515 seq = str(SeqDict[GFF[gene]['chrom']].seq[start_extension:stop_extension])
|
|
|
516
|
|
|
517 #print gene,count,pos,'\n',seq
|
|
|
518
|
|
|
519 if (seq):
|
|
|
520 res = find_stop(seq)
|
|
|
521 if res == -1 :
|
|
|
522 mylist = [gene,GFF[gene]['name'],'no stop',contexte_stop, GFF[gene]['chrom'], start_extension, stop_extension,stop_extension-start_extension,rpkm_cds,'-','-',GFF[gene]['note'],seq]
|
|
|
523 writer.writerow(mylist)
|
|
|
524 else :
|
|
|
525 indic = 'ok'
|
|
|
526 #print res
|
|
|
527 #stop_extension = start_extension + res +3
|
|
|
528 else :
|
|
10
|
529 pos = check_overlapping(gff_reader,GFF[gene]['chrom'],GFF[gene]['start']-extend,GFF[gene]['start']-1,'-',GFF[gene]['name'])
|
|
6
|
530 start_extension = pos[0]
|
|
|
531 stop_extension = pos[1]+1
|
|
10
|
532 #start_extension = GFF[gene]['start']-extend
|
|
|
533 #stop_extension = GFF[gene]['start']
|
|
6
|
534 seq = str(SeqDict[GFF[gene]['chrom']].seq[start_extension:stop_extension-1].reverse_complement())
|
|
|
535
|
|
|
536 #print gene,count,pos,'\n',seq
|
|
|
537
|
|
|
538 if (seq):
|
|
|
539 res = find_stop(seq)
|
|
|
540 if res == -1 :
|
|
|
541 mylist = [gene,GFF[gene]['name'],'no stop',contexte_stop, GFF[gene]['chrom'], start_extension, stop_extension,stop_extension-start_extension,rpkm_cds,'-','-',GFF[gene]['note'],seq]
|
|
|
542 writer.writerow(mylist)
|
|
|
543 else :
|
|
|
544 indic = 'ok'
|
|
|
545 #print res
|
|
|
546 #start_extension = stop_extension - res -3
|
|
|
547 else :
|
|
|
548 indic = 'ok'
|
|
|
549
|
|
|
550
|
|
|
551 if indic == 'ok' :
|
|
|
552 ## We save new coordinates
|
|
|
553 if GFF[gene]['strand'] == '+' :
|
|
|
554 stop_extension = start_extension + res +3
|
|
|
555 #print gene, count
|
|
|
556 #print gene,start_extension,stop_extension
|
|
|
557 seq = str(SeqDict[GFF[gene]['chrom']].seq[start_extension:stop_extension])
|
|
|
558 #print seq
|
|
|
559 count_stop = aln.count(GFF[gene]['chrom'],stop_extension-2,stop_extension+2)
|
|
|
560 if pass_length(start_extension,stop_extension) :
|
|
|
561 count_ext = aln.count(GFF[gene]['chrom'],start_extension+9,stop_extension-15)
|
|
|
562 if stop_extension > GFF[gene]['stop']+9 :
|
|
|
563 stop_ok = 1
|
|
|
564 else :
|
|
|
565 stop_ok = 0
|
|
|
566 else :
|
|
|
567 start_extension = stop_extension - res - 3
|
|
|
568 #print gene, count
|
|
|
569 #print gene,start_extension,stop_extension
|
|
|
570 seq = str(SeqDict[GFF[gene]['chrom']].seq[start_extension-1:stop_extension-1].reverse_complement())
|
|
|
571 #print seq
|
|
|
572 count_stop = aln.count(GFF[gene]['chrom'],start_extension-2,start_extension+2)
|
|
|
573 if pass_length(start_extension,stop_extension) :
|
|
|
574 count_ext = aln.count(GFF[gene]['chrom'],start_extension+15,stop_extension-9)
|
|
|
575 if start_extension < GFF[gene]['start']-9 :
|
|
|
576 stop_ok = 1
|
|
|
577 else :
|
|
|
578 stop_ok = 0
|
|
|
579 ## if we are no methionine in 5 codons following stop of CDS
|
|
|
580 if (not check_met(seq) ):
|
|
|
581 ## if we have footprint in stop codon extension and stop extension is further than cds_stop+9
|
|
|
582 if count_stop > 2 and stop_ok == 1 :
|
|
|
583 homo_cov = check_homo_coverage(gene,GFF,start_extension,stop_extension,aln)
|
|
|
584 if (homo_cov) :
|
|
|
585 '''
|
|
|
586 write result witch corresponding to readthrough
|
|
|
587 '''
|
|
|
588 ##if length of extension upper than 25 we can compute rpkm
|
|
|
589 if (not pass_length(start_extension,stop_extension)) :
|
|
|
590 len_ext = stop_extension-start_extension
|
|
|
591 rpkm_ext = 'nan'
|
|
|
592 ratio = 'nan'
|
|
|
593 else :
|
|
|
594 len_ext = stop_extension-start_extension
|
|
|
595 rpkm_ext = compute_rpkm(len_ext,count_ext,count_tot)
|
|
|
596 ## compute ratio between extension coverage and cds coverage (rpkm)
|
|
|
597 ratio = rpkm_ext/rpkm_cds
|
|
|
598 #print gene,ratio
|
|
|
599 #print start_extension,stop_extension
|
|
|
600 mylist = [gene,GFF[gene]['name'],'-',contexte_stop,GFF[gene]['chrom'], start_extension, stop_extension,stop_extension-start_extension,rpkm_cds,rpkm_ext,ratio,GFF[gene]['note'],seq]
|
|
|
601 writer.writerow(mylist)
|
|
|
602 file_context.write('>'+gene+'\n'+contexte_stop+'\n')
|
|
|
603 file_extension.write('>'+gene+'\n'+seq+'\n')
|
|
|
604 else :
|
|
|
605 '''
|
|
|
606 write result witch corresponding to readthrough but with no homogeneous coverage
|
|
|
607 '''
|
|
|
608 if (not pass_length(start_extension,stop_extension)) :
|
|
|
609 len_ext = stop_extension-start_extension
|
|
|
610 rpkm_ext = 'nan'
|
|
|
611 ratio = 'nan'
|
|
|
612 else :
|
|
|
613 len_ext = stop_extension-start_extension
|
|
|
614 rpkm_ext = compute_rpkm(len_ext,count_ext,count_tot)
|
|
|
615 ## compute ratio between extension coverage and cds coverage (rpkm)
|
|
|
616 ratio = rpkm_ext/rpkm_cds
|
|
|
617 mylist = [gene,GFF[gene]['name'],'hetero cov',contexte_stop, GFF[gene]['chrom'], start_extension, stop_extension,stop_extension-start_extension,rpkm_cds,rpkm_ext,ratio,GFF[gene]['note'],seq]
|
|
|
618 writer.writerow(mylist)
|
|
|
619 file_context.write('>'+gene+'\n'+contexte_stop+'\n')
|
|
|
620 file_extension.write('>'+gene+'\n'+seq+'\n')
|
|
|
621 #print ">"+gene+"\n"+contexte_stop
|
|
|
622
|
|
|
623 ## plot gene :
|
|
|
624 plot_gene(aln, GFF[gene], start_extension, stop_extension, dirout+"/"+gene)
|
|
|
625
|
|
|
626
|
|
|
627
|
|
|
628 else :
|
|
|
629 '''
|
|
|
630 write result with no footprint in stop codon of extension
|
|
|
631 '''
|
|
|
632 mylist = [gene,GFF[gene]['name'],'no RPF in stop',contexte_stop, GFF[gene]['chrom'], start_extension, stop_extension,stop_extension-start_extension,rpkm_cds,'-','-',GFF[gene]['note'],seq]
|
|
|
633 writer.writerow(mylist)
|
|
|
634 file_context.write('>'+gene+'\n'+contexte_stop+'\n')
|
|
|
635 file_extension.write('>'+gene+'\n'+seq+'\n')
|
|
|
636 #print ">"+gene+"\n"+contexte_stop
|
|
|
637 else :
|
|
|
638 '''
|
|
|
639 write result with RPF maybe result of reinitiation on a start codon
|
|
|
640 '''
|
|
|
641 if pass_length(start_extension,stop_extension) :
|
|
|
642 mylist = [gene,GFF[gene]['name'],'Met after stop', contexte_stop, GFF[gene]['chrom'], start_extension, stop_extension,stop_extension-start_extension,rpkm_cds,'-','-',GFF[gene]['note'],seq]
|
|
|
643 writer.writerow(mylist)
|
|
|
644 file_context.write('>'+gene+'\n'+contexte_stop+'\n')
|
|
|
645 file_extension.write('>'+gene+'\n'+seq+'\n')
|
|
|
646 #print ">"+gene+"\n"+contexte_stop
|
|
|
647 else :
|
|
|
648 ## if its not a interesting case, we get stop context of genes without readthrough
|
|
|
649 if GFF[gene]['strand'] == '+' :
|
|
|
650 contexte_stop = str(SeqDict[GFF[gene]['chrom']].seq[GFF[gene]['stop']-6:GFF[gene]['stop']+6])
|
|
|
651 file_back.write(contexte_stop+'\n')
|
|
|
652 else :
|
|
|
653 contexte_stop = str(SeqDict[GFF[gene]['chrom']].seq[GFF[gene]['start']-7:GFF[gene]['start']+5].reverse_complement())
|
|
|
654 file_back.write(contexte_stop+'\n')
|
|
|
655
|
|
|
656 ## excluded UT with incorrect positions
|
|
|
657 except ValueError:
|
|
|
658 pass
|
|
|
659
|
|
|
660
|
|
|
661 file_context.close()
|
|
|
662 file_back.close()
|
|
|
663 file_extension.close()
|
|
|
664 except Exception,e:
|
|
|
665 stop_err( 'Error during computing analysis : ' + str( e ) )
|
|
|
666
|
|
|
667
|
|
|
668 def cleaning_bam(bam):
|
|
|
669 '''
|
|
|
670 Remove reads unmapped, non uniquely mapped and reads with length lower than 25 and upper than 32, and mapping quality upper than 12
|
|
|
671 '''
|
|
|
672 try :
|
|
|
673 header = subprocess.check_output(['samtools', 'view', '-H', bam], stderr= subprocess.PIPE)
|
|
|
674 #header = results.split('\n')
|
|
|
675
|
|
|
676 # check mapper for define multiple tag
|
|
|
677 if re.search('bowtie', header):
|
|
|
678 multi_tag = "XS:i:"
|
|
|
679 elif re.search('bwa', header):
|
|
|
680 multi_tag = "XT:A:M"
|
|
|
681 else :
|
|
|
682 stop_err("No PG tag find in"+bam+". Please use bowtie or bwa for mapping")
|
|
|
683
|
|
|
684 tmp_sam = tempfile.mktemp()
|
|
|
685 cmd = "samtools view %s > %s" % (bam, tmp_sam)
|
|
|
686 proc = subprocess.Popen( args=cmd, shell=True, stderr = subprocess.PIPE)
|
|
|
687 returncode = proc.wait()
|
|
|
688
|
|
|
689
|
|
|
690 with open(tempfile.mktemp(),'w') as out :
|
|
|
691 out.write(header)
|
|
|
692 with open(tmp_sam,'r') as sam :
|
|
|
693 for line in sam :
|
|
|
694 if multi_tag not in line and line.split('\t')[1] != '4' and int(line.split('\t')[4]) > 12 :
|
|
|
695 if len(line.split('\t')[9]) < 32 and len(line.split('\t')[9]) > 25 :
|
|
|
696 out.write(line)
|
|
|
697 bamfile = tempfile.mktemp()+'.bam'
|
|
|
698 cmd = "samtools view -hSb %s > %s" % (out.name,bamfile)
|
|
|
699 proc = subprocess.Popen( args=cmd, shell=True, stderr = subprocess.PIPE)
|
|
|
700 returncode = proc.wait()
|
|
|
701
|
|
|
702 return bamfile
|
|
|
703
|
|
|
704 except Exception,e:
|
|
|
705 stop_err( 'Error during cleaning bam : ' + str( e ) )
|
|
|
706
|
|
|
707 def write_html_page(html,subfolder) :
|
|
|
708
|
|
|
709
|
|
|
710 try :
|
|
|
711
|
|
|
712 gene_table = ''
|
|
|
713 gene_table += '<table>'
|
|
|
714 gene_table += '<thead><tr><th data-sort="string">Gene</th><th>Plot</th><th data-sort="string">Name</th><th>Stop context</th><th>Coordinates</th><th>RPKM CDS</th><th>RPKM extension</th><th data-sort="float">ratio</th><th>Extension</th></tr></thead><tbody>'
|
|
|
715
|
|
|
716 with open(os.path.join(subfolder,'readthrough_result.csv'), 'rb') as csvfile:
|
|
|
717 spamreader = csv.reader(csvfile, delimiter='\t')
|
|
|
718 ## skip the header
|
|
|
719 next(spamreader, None)
|
|
|
720 for row in spamreader:
|
|
|
721 if row[2] == '-' :
|
|
|
722 gene_table += '<tr><td>%s</td><td><a href="%s.png" data-lightbox="%s"><img src="%s_thumbnail.png" /></a></td><td><a title="%s">%s</a></td><td>%s</td><td>%s:%s-%s</td><td>%s</td><td>%s</td><td>%s</td><td>%s</td></tr>' %(row[0], row[0], row[0], row[0], row[11], row[1], row[3], row[4], row[5], row[6], row[8], row[9], row[10], row[12])
|
|
|
723
|
|
|
724 gene_table += '</tbody></table>'
|
|
|
725
|
|
|
726
|
|
|
727
|
|
|
728 html_str = """
|
|
|
729 <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
|
|
|
730 "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
|
|
|
731
|
|
|
732 <html xmlns="http://www.w3.org/1999/xhtml">
|
|
|
733 <head>
|
|
|
734 <link rel="stylesheet" href="http://code.jquery.com/ui/1.10.3/themes/smoothness/jquery-ui.css" />
|
|
|
735 <script src="http://code.jquery.com/jquery-1.10.2.min.js"></script>
|
|
|
736 <script src="http://ajax.googleapis.com/ajax/libs/jquery/1.9.1/jquery.min.js"></script>
|
|
|
737 <script src="http://code.jquery.com/ui/1.10.3/jquery-ui.js"></script>
|
|
|
738 <script src="lightbox/js/jquery-1.10.2.min.js"></script>
|
|
|
739 <script src="lightbox/js/lightbox-2.6.min.js"></script>
|
|
|
740 <link href="lightbox/css/lightbox.css" rel="stylesheet" />
|
|
|
741 <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
|
|
|
742 <title>Dual coding result file</title>
|
|
|
743 <link rel="stylesheet" href="http://code.jquery.com/ui/1.10.3/themes/smoothness/jquery-ui.css" />
|
|
|
744 <script>
|
|
|
745 (function(d){d.fn.stupidtable=function(b){return this.each(function(){var a=d(this);b=b||{};b=d.extend({},d.fn.stupidtable.default_sort_fns,b);var n=function(a,b){for(var f=[],c=0,e=a.slice(0).sort(b),h=0;h<a.length;h++){for(c=d.inArray(a[h],e);-1!=d.inArray(c,f);)c++;f.push(c)}return f},q=function(a,b){for(var d=a.slice(0),c=0,e=0;e<b.length;e++)c=b[e],d[c]=a[e];return d};a.on("click","th",function(){var m=a.children("tbody").children("tr"),g=d(this),f=0,c=d.fn.stupidtable.dir;a.find("th").slice(0, g.index()).each(function(){var a=d(this).attr("colspan")||1;f+=parseInt(a,10)});var e=g.data("sort-default")||c.ASC;g.data("sort-dir")&&(e=g.data("sort-dir")===c.ASC?c.DESC:c.ASC);var h=g.data("sort")||null;null!==h&&(a.trigger("beforetablesort",{column:f,direction:e}),a.css("display"),setTimeout(function(){var l=[],p=b[h];m.each(function(a,b){var c=d(b).children().eq(f),e=c.data("sort-value"),c="undefined"!==typeof e?e:c.text();l.push(c)});var k;k=e==c.ASC?n(l,p):n(l,function(a,b){return-p(a,b)}); a.find("th").data("sort-dir",null).removeClass("sorting-desc sorting-asc");g.data("sort-dir",e).addClass("sorting-"+e);k=d(q(m,k));a.children("tbody").remove();a.append("<tbody />").append(k);a.trigger("aftertablesort",{column:f,direction:e});a.css("display")},10))})})};d.fn.stupidtable.dir={ASC:"asc",DESC:"desc"};d.fn.stupidtable.default_sort_fns={"int":function(b,a){return parseInt(b,10)-parseInt(a,10)},"float":function(b,a){return parseFloat(b)-parseFloat(a)},string:function(b,a){return b<a?-1: b>a?1:0},"string-ins":function(b,a){b=b.toLowerCase();a=a.toLowerCase();return b<a?-1:b>a?1:0}}})(jQuery);
|
|
|
746 (function($) {
|
|
|
747
|
|
|
748 $.fn.stupidtable = function(sortFns) {
|
|
|
749 return this.each(function() {
|
|
|
750 var $table = $(this);
|
|
|
751 sortFns = sortFns || {};
|
|
|
752
|
|
|
753 // ==================================================== //
|
|
|
754 // Utility functions //
|
|
|
755 // ==================================================== //
|
|
|
756
|
|
|
757 // Merge sort functions with some default sort functions.
|
|
|
758 sortFns = $.extend({}, $.fn.stupidtable.default_sort_fns, sortFns);
|
|
|
759
|
|
|
760 // Return the resulting indexes of a sort so we can apply
|
|
|
761 // this result elsewhere. This returns an array of index numbers.
|
|
|
762 // return[0] = x means "arr's 0th element is now at x"
|
|
|
763 var sort_map = function(arr, sort_function) {
|
|
|
764 var map = [];
|
|
|
765 var index = 0;
|
|
|
766 var sorted = arr.slice(0).sort(sort_function);
|
|
|
767 for (var i=0; i<arr.length; i++) {
|
|
|
768 index = $.inArray(arr[i], sorted);
|
|
|
769
|
|
|
770 // If this index is already in the map, look for the next index.
|
|
|
771 // This handles the case of duplicate entries.
|
|
|
772 while ($.inArray(index, map) != -1) {
|
|
|
773 index++;
|
|
|
774 }
|
|
|
775 map.push(index);
|
|
|
776 }
|
|
|
777
|
|
|
778 return map;
|
|
|
779 };
|
|
|
780
|
|
|
781 // Apply a sort map to the array.
|
|
|
782 var apply_sort_map = function(arr, map) {
|
|
|
783 var clone = arr.slice(0),
|
|
|
784 newIndex = 0;
|
|
|
785 for (var i=0; i<map.length; i++) {
|
|
|
786 newIndex = map[i];
|
|
|
787 clone[newIndex] = arr[i];
|
|
|
788 }
|
|
|
789 return clone;
|
|
|
790 };
|
|
|
791
|
|
|
792 // ==================================================== //
|
|
|
793 // Begin execution! //
|
|
|
794 // ==================================================== //
|
|
|
795
|
|
|
796 // Do sorting when THs are clicked
|
|
|
797 $table.on("click", "th", function() {
|
|
|
798 var trs = $table.children("tbody").children("tr");
|
|
|
799 var $this = $(this);
|
|
|
800 var th_index = 0;
|
|
|
801 var dir = $.fn.stupidtable.dir;
|
|
|
802
|
|
|
803 $table.find("th").slice(0, $this.index()).each(function() {
|
|
|
804 var cols = $(this).attr("colspan") || 1;
|
|
|
805 th_index += parseInt(cols,10);
|
|
|
806 });
|
|
|
807
|
|
|
808 // Determine (and/or reverse) sorting direction, default `asc`
|
|
|
809 var sort_dir = $this.data("sort-default") || dir.ASC;
|
|
|
810 if ($this.data("sort-dir"))
|
|
|
811 sort_dir = $this.data("sort-dir") === dir.ASC ? dir.DESC : dir.ASC;
|
|
|
812
|
|
|
813 // Choose appropriate sorting function.
|
|
|
814 var type = $this.data("sort") || null;
|
|
|
815
|
|
|
816 // Prevent sorting if no type defined
|
|
|
817 if (type === null) {
|
|
|
818 return;
|
|
|
819 }
|
|
|
820
|
|
|
821 // Trigger `beforetablesort` event that calling scripts can hook into;
|
|
|
822 // pass parameters for sorted column index and sorting direction
|
|
|
823 $table.trigger("beforetablesort", {column: th_index, direction: sort_dir});
|
|
|
824 // More reliable method of forcing a redraw
|
|
|
825 $table.css("display");
|
|
|
826
|
|
|
827 // Run sorting asynchronously on a timout to force browser redraw after
|
|
|
828 // `beforetablesort` callback. Also avoids locking up the browser too much.
|
|
|
829 setTimeout(function() {
|
|
|
830 // Gather the elements for this column
|
|
|
831 var column = [];
|
|
|
832 var sortMethod = sortFns[type];
|
|
|
833
|
|
|
834 // Push either the value of the `data-order-by` attribute if specified
|
|
|
835 // or just the text() value in this column to column[] for comparison.
|
|
|
836 trs.each(function(index,tr) {
|
|
|
837 var $e = $(tr).children().eq(th_index);
|
|
|
838 var sort_val = $e.data("sort-value");
|
|
|
839 var order_by = typeof(sort_val) !== "undefined" ? sort_val : $e.text();
|
|
|
840 column.push(order_by);
|
|
|
841 });
|
|
|
842
|
|
|
843 // Create the sort map. This column having a sort-dir implies it was
|
|
|
844 // the last column sorted. As long as no data-sort-desc is specified,
|
|
|
845 // we're free to just reverse the column.
|
|
|
846 var theMap;
|
|
|
847 if (sort_dir == dir.ASC)
|
|
|
848 theMap = sort_map(column, sortMethod);
|
|
|
849 else
|
|
|
850 theMap = sort_map(column, function(a, b) { return -sortMethod(a, b); });
|
|
|
851
|
|
|
852 // Reset siblings
|
|
|
853 $table.find("th").data("sort-dir", null).removeClass("sorting-desc sorting-asc");
|
|
|
854 $this.data("sort-dir", sort_dir).addClass("sorting-"+sort_dir);
|
|
|
855
|
|
|
856 var sortedTRs = $(apply_sort_map(trs, theMap));
|
|
|
857 $table.children("tbody").remove();
|
|
|
858 $table.append("<tbody />").append(sortedTRs);
|
|
|
859
|
|
|
860 // Trigger `aftertablesort` event. Similar to `beforetablesort`
|
|
|
861 $table.trigger("aftertablesort", {column: th_index, direction: sort_dir});
|
|
|
862 // More reliable method of forcing a redraw
|
|
|
863 $table.css("display");
|
|
|
864 }, 10);
|
|
|
865 });
|
|
|
866 });
|
|
|
867 };
|
|
|
868
|
|
|
869 // Enum containing sorting directions
|
|
|
870 $.fn.stupidtable.dir = {ASC: "asc", DESC: "desc"};
|
|
|
871
|
|
|
872 $.fn.stupidtable.default_sort_fns = {
|
|
|
873 "int": function(a, b) {
|
|
|
874 return parseInt(a, 10) - parseInt(b, 10);
|
|
|
875 },
|
|
|
876 "float": function(a, b) {
|
|
|
877 return parseFloat(a) - parseFloat(b);
|
|
|
878 },
|
|
|
879 "string": function(a, b) {
|
|
|
880 if (a < b) return -1;
|
|
|
881 if (a > b) return +1;
|
|
|
882 return 0;
|
|
|
883 },
|
|
|
884 "string-ins": function(a, b) {
|
|
|
885 a = a.toLowerCase();
|
|
|
886 b = b.toLowerCase();
|
|
|
887 if (a < b) return -1;
|
|
|
888 if (a > b) return +1;
|
|
|
889 return 0;
|
|
|
890 }
|
|
|
891 };
|
|
|
892
|
|
|
893 })(jQuery);
|
|
|
894 $(function(){
|
|
|
895 var table = $("table").stupidtable();
|
|
|
896
|
|
|
897 table.on("beforetablesort", function (event, data) {
|
|
|
898 // data.column - the index of the column sorted after a click
|
|
|
899 // data.direction - the sorting direction (either asc or desc)
|
|
|
900 $("#msg").text("Sorting index " + data.column)
|
|
|
901 });
|
|
|
902 table.on("aftertablesort", function (event, data) {
|
|
|
903 var th = $(this).find("th");
|
|
|
904 th.find(".arrow").remove();
|
|
|
905 var dir = $.fn.stupidtable.dir;
|
|
|
906 var arrow = data.direction === dir.ASC ? "↑" : "↓";
|
|
|
907 th.eq(data.column).append('<span class="arrow">' + arrow +'</span>');
|
|
|
908 });
|
|
|
909 });
|
|
|
910 </script>
|
|
|
911 <style type="text/css">
|
|
|
912 label {
|
|
|
913 display: inline-block;
|
|
|
914 width: 5em;
|
|
|
915 }
|
|
|
916 table {
|
|
|
917 border-collapse: collapse;
|
|
|
918 }
|
|
|
919 th, td {
|
|
|
920 padding: 5px 10px;
|
|
|
921 border: 1px solid #999;
|
|
|
922 }
|
|
|
923 th {
|
|
|
924 background-color: #a7d3ff;
|
|
|
925 }
|
|
|
926 th[data-sort]{
|
|
|
927 cursor:pointer;
|
|
|
928 }
|
|
|
929 a[data-lightbox]{
|
|
|
930 cursor:zoom-in;
|
|
|
931 }
|
|
|
932 #msg {
|
|
|
933 color: green;
|
|
|
934 }
|
|
|
935 </style>
|
|
|
936 </head>
|
|
|
937 <body>
|
|
|
938 <h1>Readthrough analyse results</h1>
|
|
|
939 %s
|
|
|
940 </body>
|
|
|
941 </html> """ % (gene_table)
|
|
|
942
|
|
|
943 html_file = open(html,"w")
|
|
|
944 html_file.write(html_str)
|
|
|
945 html_file.close()
|
|
|
946
|
|
|
947
|
|
|
948 except Exception, e :
|
|
|
949 stop_err('Error during html page creation : ' + str( e ) )
|
|
|
950
|
|
|
951
|
|
|
952 def __main__():
|
|
|
953
|
|
|
954 #Parse command line options
|
|
|
955 parser = optparse.OptionParser()
|
|
|
956 parser.add_option("-g", "--gff", dest="gff", type= "string",
|
|
|
957 help="GFF annotation file", metavar="FILE")
|
|
|
958
|
|
|
959 parser.add_option("-f", "--fasta", dest="fasta", type= "string",
|
|
|
960 help="Fasta file ", metavar="FILE")
|
|
|
961
|
|
|
962 parser.add_option("-b", "--bam", dest="bamfile", type= "string",
|
|
|
963 help="Bam Ribo-Seq alignments ", metavar="FILE")
|
|
|
964
|
|
|
965 parser.add_option("-d", "--dirout", dest="dirout", type= "string",
|
|
10
|
966 help="write report in this html file and in associated directory", metavar="STR,STR")
|
|
|
967
|
|
|
968 parser.add_option("-e", "--extend", dest="extend", type= "int",
|
|
|
969 help="Lenght of extension after stop in number of base pairs (depends on your organisme)", metavar="INT")
|
|
6
|
970
|
|
|
971 parser.add_option("-q", "--quiet",
|
|
|
972 action="store_false", dest="verbose", default=True,
|
|
|
973 help="don't print status messages to stdout")
|
|
|
974
|
|
|
975 (options, args) = parser.parse_args()
|
|
|
976 sys.stdout.write("Begin readthrough analysis at %s\n" % time.asctime( time.localtime(time.time())))
|
|
|
977
|
|
|
978 try:
|
|
|
979 (html_file, subfolder ) = options.dirout.split(",")
|
|
|
980 if os.path.exists(subfolder):
|
|
|
981 raise
|
|
|
982 try:
|
|
|
983 os.mkdir(subfolder)
|
|
|
984 except:
|
|
|
985 raise
|
|
10
|
986 ## identify GFF or GTF format from 9th column
|
|
|
987 with open (options.gff,"r") as gffile :
|
|
|
988 for line in gffile :
|
|
|
989 if '#' in line :
|
|
|
990 ## skip header
|
|
|
991 gffile.next()
|
|
|
992 elif 'gene_id' in line :
|
|
|
993 ## launch gtf reader :
|
|
|
994 GFF = ribo_functions.store_gtf(options.gff)
|
|
|
995 break
|
|
|
996 elif 'ID=' in line :
|
|
|
997 ## launch gff reader
|
|
|
998 GFF = ribo_functions.store_gff(options.gff)
|
|
|
999 break
|
|
|
1000 else :
|
|
|
1001 stop_err( 'Please check your annotation file is in correct format, GFF or GTF' )
|
|
|
1002
|
|
|
1003 #GFF = store_gff(options.gff)
|
|
|
1004 #GFF = ribo_functions.store_gtf(options.gff)
|
|
|
1005 ## check gff reading
|
|
|
1006 if not GFF['order'] :
|
|
|
1007 stop_err( 'Incorrect GFF file' + str( e ) )
|
|
6
|
1008 clean_file = cleaning_bam(options.bamfile)
|
|
10
|
1009 compute_analysis(clean_file, GFF, options.fasta, options.gff, subfolder, options.extend)
|
|
6
|
1010 if os.path.exists( clean_file ):
|
|
|
1011 os.remove( clean_file )
|
|
|
1012
|
|
|
1013 write_html_page(html_file,subfolder)
|
|
|
1014 ##paste jquery script in result directory :
|
|
|
1015 jq_src = os.path.join(os.path.dirname(__file__),'lightbox')
|
|
|
1016 shutil.copytree(jq_src,os.path.join(subfolder,'lightbox'))
|
|
|
1017
|
|
|
1018
|
|
|
1019 sys.stdout.write("Finish readthrough analysis at %s\n" % time.asctime( time.localtime(time.time())))
|
|
|
1020 except Exception, e:
|
|
|
1021 stop_err( 'Error running metagene readthrough analysis : ' + str( e ) )
|
|
|
1022
|
|
|
1023
|
|
|
1024 if __name__=="__main__":
|
|
|
1025 __main__()
|
