|
9
|
1 #!/usr/bin/env python2.7
|
|
|
2
|
|
|
3 '''
|
|
|
4 Created on Jan. 2014
|
|
|
5 @author: rachel legendre
|
|
|
6 @copyright: rachel.legendre@igmors.u-psud.fr
|
|
|
7 @license: GPL v3
|
|
|
8 '''
|
|
|
9
|
|
13
|
10 import sys, subprocess, re, commands, time, urllib, tempfile
|
|
10
|
11 from copy import copy
|
|
|
12
|
|
9
|
13
|
|
|
14 def stop_err( msg ):
|
|
|
15 sys.stderr.write( "%s\n" % msg )
|
|
|
16 sys.stderr.write( "Programme aborted at %s\n" % time.asctime(time.localtime(time.time())))
|
|
|
17 sys.exit()
|
|
|
18
|
|
|
19 def split_bam(bamfile,tmpdir):
|
|
|
20 '''
|
|
|
21 split bam by chromosome and write sam file in tmpdir
|
|
|
22 '''
|
|
|
23 try:
|
|
|
24 #get header
|
|
|
25 results = subprocess.check_output(['samtools', 'view', '-H',bamfile])
|
|
|
26 header = results.split('\n')
|
|
|
27
|
|
|
28 #define genome size
|
|
|
29 genome = []
|
|
|
30 for line in header:
|
|
|
31 result = re.search('SN', line)
|
|
|
32 if result :
|
|
|
33 #print line
|
|
|
34 feat = line.split('\t')
|
|
|
35 chrom = re.split(":", feat[1])
|
|
|
36 #print feat[1]
|
|
|
37 genome.append(chrom[1])
|
|
|
38
|
|
|
39 #split sam by chrom
|
|
|
40 n = 0
|
|
|
41 for chrm in genome:
|
|
|
42 with open(tmpdir+'/'+chrm+'.sam', 'w') as f :
|
|
|
43 #write header correctly for each chromosome
|
|
|
44 f.write(header[0]+'\n')
|
|
|
45 expr = re.compile(chrm+'\t')
|
|
|
46 el =[elem for elem in header if expr.search(elem)][0]
|
|
|
47 f.write(el+'\n')
|
|
|
48 f.write(header[-2]+'\n')
|
|
|
49 #write all reads for each chromosome
|
|
|
50 reads = subprocess.check_output(["samtools", "view", bamfile, chrm])
|
|
|
51 f.write(reads)
|
|
|
52 # calculate number of reads
|
|
|
53 n += reads.count(chrm)
|
|
|
54
|
|
|
55 sys.stdout.write("%d reads are presents in your bam file\n" % n)
|
|
|
56
|
|
|
57 except Exception, e:
|
|
|
58 stop_err( 'Error during bam file splitting : ' + str( e ) )
|
|
|
59
|
|
|
60
|
|
|
61
|
|
17
|
62 def get_first_base(tmpdir, kmer, frame):
|
|
9
|
63 '''
|
|
|
64 write footprint coverage file for each sam file in tmpdir
|
|
|
65 '''
|
|
|
66 global total_mapped_read
|
|
17
|
67 ## tags by default
|
|
|
68 multi_tag = "XS:i:"
|
|
|
69 tag = "IH:i:1"
|
|
|
70
|
|
9
|
71 try:
|
|
17
|
72 ### compute position of P-site according to frame (P-site -> +15 starting from 5' end of footprint)
|
|
|
73 p_site_pos = 16-frame
|
|
|
74
|
|
9
|
75 file_array = (commands.getoutput('ls '+tmpdir)).split('\n')
|
|
|
76 ##write coverage for each sam file in tmpdir
|
|
|
77 for samfile in file_array:
|
|
|
78 with open(tmpdir+'/'+samfile, 'r') as sam :
|
|
|
79 ##get chromosome name
|
|
17
|
80 chrom = samfile.split(".sam")[0]
|
|
9
|
81
|
|
|
82 for line in sam:
|
|
|
83 #initialize dictionnary
|
|
17
|
84 if '@SQ' in line :
|
|
9
|
85 size = int(line.split('LN:')[1])
|
|
|
86 genomeF = [0]*size
|
|
|
87 genomeR = [0]*size
|
|
|
88 # define multiple reads keys from mapper
|
|
17
|
89 elif '@PG\tID' in line :
|
|
|
90 if 'bowtie' in line:
|
|
9
|
91 multi_tag = "XS:i:"
|
|
17
|
92 elif 'bwa' in line:
|
|
9
|
93 multi_tag = "XT:A:R"
|
|
17
|
94 elif 'TopHat' in line:
|
|
|
95 tag = "NH:i:1"
|
|
9
|
96 else :
|
|
17
|
97 stop_err("No PG tag find in "+samfile+". Please use bowtie, bwa or Tophat for mapping")
|
|
9
|
98
|
|
|
99 # get footprint
|
|
|
100 elif re.search('^[^@].+', line) :
|
|
|
101 #print line.rstrip()
|
|
|
102 read_pos = int(line.split('\t')[3])
|
|
|
103 read_sens = int(line.split('\t')[1])
|
|
17
|
104 len_read = len(line.split('\t')[9])
|
|
|
105 read_qual = int(line.split('\t')[4])
|
|
|
106 if len_read == kmer and (tag in line or multi_tag not in line) and read_qual > 20 :
|
|
9
|
107 ###if line.split('\t')[5] == '28M' :
|
|
|
108 # print line.rstrip()
|
|
|
109 total_mapped_read +=1
|
|
|
110 #if it's a forward read
|
|
|
111 if read_sens == 0 :
|
|
|
112 #get P-site : start read + 12 nt
|
|
17
|
113 #read_pos += 15
|
|
|
114 read_pos += p_site_pos
|
|
9
|
115 genomeF[read_pos] += 1
|
|
|
116 #if it's a reverse read
|
|
|
117 elif read_sens == 16 :
|
|
|
118 #get P-site
|
|
17
|
119 #read_pos += 12
|
|
|
120 read_pos += (len_read-1) - p_site_pos
|
|
9
|
121 genomeR[read_pos] += 1
|
|
|
122
|
|
|
123 #try:
|
|
|
124 #write coverage in files
|
|
|
125 with open(tmpdir+'/assoCov_'+chrom+'.txt', 'w') as cov :
|
|
|
126 for i in range(0,size):
|
|
|
127 cov.write(str(genomeF[i])+'\t-'+str(genomeR[i])+'\n')
|
|
|
128 #except Exception,e:
|
|
|
129 # stop_err( 'Error during coverage file writting : ' + str( e ) )
|
|
17
|
130 del genomeF
|
|
|
131 del genomeR
|
|
9
|
132 sys.stdout.write("%d reads are used for frame analysis\n" % total_mapped_read)
|
|
|
133 except Exception, e:
|
|
|
134 stop_err( 'Error during footprinting : ' + str( e ) )
|
|
17
|
135
|
|
|
136
|
|
9
|
137
|
|
|
138 def store_gff(gff):
|
|
|
139 '''
|
|
|
140 parse and store gff file in a dictionnary
|
|
|
141 '''
|
|
|
142 try:
|
|
|
143 GFF = {}
|
|
18
|
144 mRNA = {}
|
|
9
|
145 with open(gff, 'r') as f_gff :
|
|
13
|
146
|
|
9
|
147 GFF['order'] = []
|
|
13
|
148 #for line in itertools.islice( f_gff, 25 ):
|
|
9
|
149 for line in f_gff:
|
|
|
150 ## switch commented lines
|
|
|
151 line = line.split("#")[0]
|
|
|
152 if line != "" :
|
|
|
153 feature = (line.split('\t')[8]).split(';')
|
|
|
154 # first line is already gene line :
|
|
|
155 if line.split('\t')[2] == 'gene' :
|
|
|
156 gene = feature[0].replace("ID=","")
|
|
18
|
157 curent_gene = gene
|
|
13
|
158 if 'Name' in line :
|
|
18
|
159 regex = re.compile('(Name=)([^;]*);?')
|
|
13
|
160 res = regex.search(line.split('\t')[8])
|
|
|
161 Name = res.group(2)
|
|
|
162 Name = Name.rstrip()
|
|
9
|
163 else :
|
|
|
164 Name = "Unknown"
|
|
|
165 ##get annotation
|
|
13
|
166 if 'Note' in line :
|
|
18
|
167 regex = re.compile('(Note=)([^;]*);?')
|
|
13
|
168 res = regex.search(line.split('\t')[8])
|
|
|
169 note = res.group(2)
|
|
|
170 note = urllib.unquote(str(note)).replace("\n","")
|
|
|
171 else :
|
|
|
172 note = ""
|
|
9
|
173 ## store gene information
|
|
|
174 GFF['order'].append(gene)
|
|
|
175 GFF[gene] = {}
|
|
|
176 GFF[gene]['chrom'] = line.split('\t')[0]
|
|
|
177 GFF[gene]['start'] = int(line.split('\t')[3])
|
|
|
178 GFF[gene]['stop'] = int(line.split('\t')[4])
|
|
|
179 GFF[gene]['strand'] = line.split('\t')[6]
|
|
|
180 GFF[gene]['name'] = Name
|
|
|
181 GFF[gene]['note'] = note
|
|
|
182 GFF[gene]['exon'] = {}
|
|
|
183 GFF[gene]['exon_number'] = 0
|
|
|
184 #print Name
|
|
18
|
185 elif line.split('\t')[2] == 'mRNA' :
|
|
|
186 regex = re.compile('(Parent=)([^;]*);?')
|
|
|
187 res = regex.search(line.split('\t')[8])
|
|
|
188 gene_name = res.group(2)
|
|
|
189 regex = re.compile('(ID=)([^;]*);?')
|
|
|
190 res = regex.search(line.split('\t')[8])
|
|
|
191 mRNA_name = res.group(2)
|
|
|
192 if gene not in mRNA.viewvalues() and gene_name == curent_gene :
|
|
|
193 mRNA[mRNA_name] = gene_name
|
|
|
194
|
|
9
|
195 elif line.split('\t')[2] == 'CDS' :
|
|
18
|
196 regex = re.compile('(Parent=)([^;]*);?')
|
|
13
|
197 res = regex.search(line.split('\t')[8])
|
|
|
198 gene = res.group(2)
|
|
18
|
199
|
|
13
|
200 if 'mRNA' in gene:
|
|
|
201 gene = re.sub(r"(.*)(\_mRNA)", r"\1", gene)
|
|
18
|
202 if mRNA.has_key(gene) and GFF.has_key(mRNA[gene]):
|
|
|
203
|
|
|
204 gene = gene_name
|
|
|
205
|
|
9
|
206 if GFF.has_key(gene) :
|
|
|
207 GFF[gene]['exon_number'] += 1
|
|
|
208 exon_number = GFF[gene]['exon_number']
|
|
|
209 GFF[gene]['exon'][exon_number] = {}
|
|
|
210 GFF[gene]['exon'][exon_number]['frame'] = line.split('\t')[7]
|
|
|
211 GFF[gene]['exon'][exon_number]['start'] = int(line.split('\t')[3])
|
|
|
212 GFF[gene]['exon'][exon_number]['stop'] = int(line.split('\t')[4])
|
|
18
|
213
|
|
9
|
214 ## if there is a five prim UTR intron, we change start of gene
|
|
13
|
215 elif line.split('\t')[2] == 'five_prime_UTR_intron' :
|
|
9
|
216 if GFF[gene]['strand'] == "+" :
|
|
|
217 GFF[gene]['start'] = GFF[gene]['exon'][1]['start']
|
|
|
218 else :
|
|
|
219 GFF[gene]['stop'] = GFF[gene]['exon'][exon_number]['stop']
|
|
|
220 return GFF
|
|
|
221 except Exception,e:
|
|
|
222 stop_err( 'Error during gff storage : ' + str( e ) )
|
|
|
223
|
|
|
224
|
|
|
225 #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
|
|
|
226 #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
|
|
|
227 #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
|
|
|
228 #om%20the%20whole%20genome%20duplication;display=Vesicle%20membrane%20receptor%20protein%20%28v-SNARE%29;dbxref=SGD:S000000028;orf_classification=Verified
|
|
|
229 #chrI SGD CDS 87286 87387 . + 0 Parent=YAL030W_mRNA;Name=YAL030W_CDS;orf_classification=Verified
|
|
|
230 #chrI SGD CDS 87501 87752 . + 0 Parent=YAL030W_mRNA;Name=YAL030W_CDS;orf_classification=Verified
|
|
|
231 #chrI SGD intron 87388 87500 . + . Parent=YAL030W_mRNA;Name=YAL030W_intron;orf_classification=Verified
|
|
|
232
|
|
|
233 def store_gtf(gff):
|
|
|
234 '''
|
|
10
|
235 parse and store gtf file in a dictionnary
|
|
9
|
236 '''
|
|
|
237 try:
|
|
|
238 GFF = {}
|
|
|
239 with open(gff, 'r') as f_gff :
|
|
|
240 GFF['order'] = []
|
|
|
241 for line in f_gff:
|
|
|
242 ## switch commented lines
|
|
|
243 line = line.split("#")[0]
|
|
|
244 if line != "" :
|
|
|
245 # first line is already gene line :
|
|
10
|
246 if 'protein_coding' in line :
|
|
9
|
247 ##get name
|
|
13
|
248 gene = re.sub(r".+transcript_id \"([\w|-|\.]+)\";.*", r"\1", line).rstrip()
|
|
10
|
249 Name = re.sub(r".+gene_name \"([\w|\-|\:|\.|\(|\)]+)\";.*", r"\1", line).rstrip()
|
|
|
250 if line.split('\t')[2] == 'CDS' :
|
|
9
|
251 ##if its first time we get this gene
|
|
|
252 if gene not in GFF.keys() :
|
|
|
253 ## store gene information
|
|
|
254 GFF['order'].append(gene)
|
|
|
255 GFF[gene] = {}
|
|
|
256 GFF[gene]['chrom'] = line.split('\t')[0]
|
|
|
257 GFF[gene]['strand'] = line.split('\t')[6]
|
|
10
|
258 GFF[gene]['start'] = int(line.split('\t')[3])
|
|
|
259 GFF[gene]['stop'] = int(line.split('\t')[4])
|
|
9
|
260 GFF[gene]['name'] = Name
|
|
10
|
261 GFF[gene]['note'] = ""
|
|
9
|
262 GFF[gene]['exon_number'] = 1
|
|
|
263 GFF[gene]['exon'] = {}
|
|
10
|
264 #exon_number = int(re.sub(r".+exon_number \"(\d+)\".+", r"\1",line).rstrip())
|
|
|
265 ## some exons are non codant
|
|
|
266 exon_number = 1
|
|
9
|
267 GFF[gene]['exon'][exon_number] = {}
|
|
|
268 GFF[gene]['exon'][exon_number]['start'] = int(line.split('\t')[3])
|
|
|
269 GFF[gene]['exon'][exon_number]['stop'] = int(line.split('\t')[4])
|
|
|
270 else :
|
|
|
271 ## we add exon
|
|
10
|
272 #exon_number = int(re.sub(r".+exon_number \"(\d+)\".+", r"\1",line).rstrip())
|
|
|
273 exon_number += 1
|
|
9
|
274 GFF[gene]['exon_number'] = exon_number
|
|
|
275 GFF[gene]['exon'][exon_number] = {}
|
|
|
276 GFF[gene]['exon'][exon_number]['start'] = int(line.split('\t')[3])
|
|
|
277 GFF[gene]['exon'][exon_number]['stop'] = int(line.split('\t')[4])
|
|
10
|
278 #elif line.split('\t')[2] == 'CDS' :
|
|
|
279 #exon_number = int(re.sub(r".+exon_number \"(\d+)\".+", r"\1",line).rstrip())
|
|
9
|
280 GFF[gene]['exon'][exon_number]['frame'] = line.split('\t')[7]
|
|
|
281 elif line.split('\t')[2] == 'start_codon' :
|
|
|
282 if GFF[gene]['strand'] == '-' :
|
|
10
|
283 GFF[gene]['stop'] = int(line.split('\t')[4])
|
|
9
|
284 else :
|
|
|
285 GFF[gene]['start'] = int(line.split('\t')[3])
|
|
|
286 elif line.split('\t')[2] == 'stop_codon' :
|
|
|
287 if GFF[gene]['strand'] == '-' :
|
|
10
|
288 GFF[gene]['start'] = int(line.split('\t')[3])
|
|
9
|
289 else :
|
|
|
290 GFF[gene]['stop'] = int(line.split('\t')[4])
|
|
|
291
|
|
10
|
292 return __reverse_coordinates__(GFF)
|
|
9
|
293 except Exception,e:
|
|
13
|
294 stop_err( 'Error during gtf storage : ' + str( e ) )
|
|
9
|
295
|
|
|
296
|
|
|
297 ##IV protein_coding exon 307766 307789 . - . gene_id "YDL083C"; transcript_id "YDL083C"; exon_number "1"; gene_name "RPS16B"; gene_biotype "protein_coding"; transcript_name "RPS16B";
|
|
|
298 ## exon_id "YDL083C.1";
|
|
|
299 ##IV protein_coding CDS 307766 307789 . - 0 gene_id "YDL083C"; transcript_id "YDL083C"; exon_number "1"; gene_name "RPS16B"; gene_biotype "protein_coding"; transcript_name "RPS16B";
|
|
|
300 ## protein_id "YDL083C";
|
|
|
301 ##IV protein_coding start_codon 307787 307789 . - 0 gene_id "YDL083C"; transcript_id "YDL083C"; exon_number "1"; gene_name "RPS16B"; gene_biotype "protein_coding"; transcript_name "
|
|
|
302 ##RPS16B";
|
|
|
303 ##IV protein_coding exon 306926 307333 . - . gene_id "YDL083C"; transcript_id "YDL083C"; exon_number "2"; gene_name "RPS16B"; gene_biotype "protein_coding"; transcript_name "RPS16B";
|
|
|
304 ## exon_id "YDL083C.2";
|
|
|
305 ##IV protein_coding CDS 306929 307333 . - 0 gene_id "YDL083C"; transcript_id "YDL083C"; exon_number "2"; gene_name "RPS16B"; gene_biotype "protein_coding"; transcript_name "RPS16B";
|
|
|
306 ## protein_id "YDL083C";
|
|
|
307 ##IV protein_coding stop_codon 306926 306928 . - 0 gene_id "YDL083C"; transcript_id "YDL083C"; exon_number "2"; gene_name "RPS16B"; gene_biotype "protein_coding"; transcript_name "
|
|
|
308 ##RPS16B";
|
|
10
|
309 def __reverse_coordinates__(GFF):
|
|
|
310
|
|
|
311 for gene in GFF['order']:
|
|
|
312 ## for reverse gene
|
|
|
313 if GFF[gene]['strand'] == "-":
|
|
|
314 ## if this gene have many exon and the stop of gene is the stop of first (and not last) exon, we reverse exon coordinates
|
|
|
315 if GFF[gene]['stop'] == GFF[gene]['exon'][1]['stop'] and GFF[gene]['exon_number'] > 1 :
|
|
|
316 tmp = copy(GFF[gene]['exon'])
|
|
|
317 exon_number = GFF[gene]['exon_number']
|
|
|
318 rev_index = exon_number+1
|
|
|
319 for z in range(1,exon_number+1):
|
|
|
320 rev_index -= 1
|
|
|
321 GFF[gene]['exon'][z] = tmp[rev_index]
|
|
9
|
322
|
|
10
|
323 ## check start
|
|
|
324 if GFF[gene]['start'] != GFF[gene]['exon'][1]['start'] and GFF[gene]['start']:
|
|
|
325 GFF[gene]['exon'][1]['start'] = GFF[gene]['start']
|
|
9
|
326
|
|
10
|
327 return GFF
|
|
|
328
|
|
9
|
329
|
|
|
330 def cleaning_bam(bam):
|
|
|
331 '''
|
|
|
332 Remove reads unmapped, non uniquely mapped and reads with length lower than 25 and upper than 32, and mapping quality upper than 12
|
|
|
333 '''
|
|
|
334 try :
|
|
|
335 header = subprocess.check_output(['samtools', 'view', '-H', bam], stderr= subprocess.PIPE)
|
|
|
336 #header = results.split('\n')
|
|
13
|
337 ## tags by default
|
|
|
338 multi_tag = "XS:i:"
|
|
|
339 tag = "IH:i:1"
|
|
9
|
340 # check mapper for define multiple tag
|
|
13
|
341 if 'bowtie' in header:
|
|
9
|
342 multi_tag = "XS:i:"
|
|
13
|
343 elif 'bwa' in header:
|
|
9
|
344 multi_tag = "XT:A:R"
|
|
13
|
345 elif 'TopHat' in header:
|
|
|
346 tag = "NH:i:1"
|
|
9
|
347 else :
|
|
13
|
348 stop_err("No PG tag find in "+samfile+". Please use bowtie, bwa or Tophat for mapping")
|
|
|
349
|
|
9
|
350 tmp_sam = tempfile.mktemp()
|
|
|
351 cmd = "samtools view %s > %s" % (bam, tmp_sam)
|
|
|
352 proc = subprocess.Popen( args=cmd, shell=True, stderr = subprocess.PIPE)
|
|
|
353 returncode = proc.wait()
|
|
|
354
|
|
|
355
|
|
|
356 with open(tempfile.mktemp(),'w') as out :
|
|
|
357 out.write(header)
|
|
|
358 with open(tmp_sam,'r') as sam :
|
|
|
359 for line in sam :
|
|
13
|
360 if (multi_tag not in line or tag in line) and line.split('\t')[1] != '4' and int(line.split('\t')[4]) > 12 :
|
|
9
|
361 if len(line.split('\t')[9]) < 32 and len(line.split('\t')[9]) > 25 :
|
|
|
362 out.write(line)
|
|
|
363 bamfile = tempfile.mktemp()+'.bam'
|
|
|
364 cmd = "samtools view -hSb %s > %s" % (out.name,bamfile)
|
|
|
365 proc = subprocess.Popen( args=cmd, shell=True, stderr = subprocess.PIPE)
|
|
|
366 returncode = proc.wait()
|
|
|
367
|
|
|
368 return bamfile
|
|
|
369
|
|
|
370 except Exception,e:
|
|
|
371 stop_err( 'Error during cleaning bam : ' + str( e ) )
|
|
|
372 |
