Mercurial > repos > jackcurragh > trips_viz_bam_to_sqlite

comparison trips_bam_to_sqlite_builtin/bam_to_sqlite.py @ 13:df4bb52b226d draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author jackcurragh
date Thu, 03 Nov 2022 12:25:58 +0000
parents
children
comparison
equal deleted inserted replaced
12:02874b1b2015 13:df4bb52b226d
1 import sys
2 import pysam
3 import operator
4 import os
5 import time
6 import sqlite3
7 from sqlitedict import SqliteDict
8
9 def tran_to_genome(tran, pos, transcriptome_info_dict):
10 #print ("tran",list(transcriptome_info_dict))
11 traninfo = transcriptome_info_dict[tran]
12 chrom = traninfo["chrom"]
13 strand = traninfo["strand"]
14 exons = sorted(traninfo["exons"])
15 #print exons
16 if strand == "+":
17 exon_start = 0
18 for tup in exons:
19 exon_start = tup[0]
20 exonlen = tup[1] - tup[0]
21 if pos > exonlen:
22 pos = (pos - exonlen)-1
23 else:
24 break
25 genomic_pos = (exon_start+pos)-1
26 elif strand == "-":
27 exon_start = 0
28 for tup in exons[::-1]:
29 exon_start = tup[1]
30 exonlen = tup[1] - tup[0]
31 if pos > exonlen:
32 pos = (pos - exonlen)-1
33 else:
34 break
35 genomic_pos = (exon_start-pos)+1
36 return (chrom, genomic_pos)
37
38
39 # Takes a dictionary with a readname as key and a list of lists as value, each sub list has consists of two elements a transcript and the position the read aligns to in the transcript
40 # This function will count the number of genes that the transcripts correspond to and if less than or equal to 3 will add the relevant value to transcript_counts_dict
41 def processor(process_chunk, master_read_dict, transcriptome_info_dict,master_dict,readseq, unambig_read_length_dict):
42 readlen = len(readseq)
43 ambiguously_mapped_reads = 0
44 #get the read name
45 read = list(process_chunk)[0]
46
47 read_list = process_chunk[read] # a list of lists of all transcripts the read aligns to and the positions
48 #used to store different genomic poistions
49 genomic_positions = []
50
51 #This section is just to get the different genomic positions the read aligns to
52
53 for listname in process_chunk[read]:
54
55 tran = listname[0].replace("-","_").replace("(","").replace(")","")
56
57 pos = int(listname[1])
58 genomic_pos = tran_to_genome(tran, pos, transcriptome_info_dict)
59 #print ("genomic pos",genomic_pos)
60 if genomic_pos not in genomic_positions:
61 genomic_positions.append(genomic_pos)
62
63 #If the read maps unambiguously
64 if len(genomic_positions) == 1:
65 if readlen not in unambig_read_length_dict:
66 unambig_read_length_dict[readlen] = 0
67 unambig_read_length_dict[readlen] += 1
68 #assume this read aligns to a noncoding position, if we find that it does align to a coding region change this to True
69 coding=False
70
71 # For each transcript this read alings to
72 for listname in process_chunk[read]:
73 #get the transcript name
74 tran = listname[0].replace("-","_").replace("(","").replace(")","")
75 #If we haven't come across this transcript already then add to master_read_dict
76 if tran not in master_read_dict:
77 master_read_dict[tran] = {"ambig":{}, "unambig":{}, "mismatches":{}, "seq":{}}
78 #get the raw unedited positon, and read tags
79 pos = int(listname[1])
80 read_tags = listname[2]
81 #If there is mismatches in this line, then modify the postion and readlen (if mismatches at start or end) and add mismatches to dictionary
82 nm_tag = 0
83
84 for tag in read_tags:
85 if tag[0] == "NM":
86 nm_tag = int(tag[1])
87 if nm_tag > 0:
88 md_tag = ""
89 for tag in read_tags:
90 if tag[0] == "MD":
91 md_tag = tag[1]
92 pos_modifier, readlen_modifier,mismatches = get_mismatch_pos(md_tag,pos,readlen,master_read_dict,tran,readseq)
93 # Count the mismatches (we only do this for unambiguous)
94 for mismatch in mismatches:
95 #Ignore mismatches appearing in the first position (due to non templated addition)
96 if mismatch != 0:
97 char = mismatches[mismatch]
98 mismatch_pos = pos + mismatch
99 if mismatch_pos not in master_read_dict[tran]["seq"]:
100 master_read_dict[tran]["seq"][mismatch_pos] = {}
101 if char not in master_read_dict[tran]["seq"][mismatch_pos]:
102 master_read_dict[tran]["seq"][mismatch_pos][char] = 0
103 master_read_dict[tran]["seq"][mismatch_pos][char] += 1
104 # apply the modifiers
105 #pos = pos+pos_modifier
106 #readlen = readlen - readlen_modifier
107
108
109 try:
110 cds_start = transcriptome_info_dict[tran]["cds_start"]
111 cds_stop = transcriptome_info_dict[tran]["cds_stop"]
112
113 if pos >= cds_start and pos <= cds_stop:
114 coding=True
115 except:
116 pass
117
118
119 if readlen in master_read_dict[tran]["unambig"]:
120 if pos in master_read_dict[tran]["unambig"][readlen]:
121 master_read_dict[tran]["unambig"][readlen][pos] += 1
122 else:
123 master_read_dict[tran]["unambig"][readlen][pos] = 1
124 else:
125 master_read_dict[tran]["unambig"][readlen] = {pos:1}
126
127 if coding == True:
128 master_dict["unambiguous_coding_count"] += 1
129 elif coding == False:
130 master_dict["unambiguous_non_coding_count"] += 1
131
132 else:
133 ambiguously_mapped_reads += 1
134 for listname in process_chunk[read]:
135 tran = listname[0].replace("-","_").replace("(","").replace(")","")
136 if tran not in master_read_dict:
137 master_read_dict[tran] = {"ambig":{}, "unambig":{}, "mismatches":{}, "seq":{}}
138 pos = int(listname[1])
139 read_tags = listname[2]
140 nm_tag = 0
141 for tag in read_tags:
142 if tag[0] == "NM":
143 nm_tag = int(tag[1])
144 if nm_tag > 0:
145 md_tag = ""
146 for tag in read_tags:
147 if tag[0] == "MD":
148 md_tag = tag[1]
149 pos_modifier, readlen_modifier,mismatches = get_mismatch_pos(md_tag,pos,readlen,master_read_dict,tran,readseq)
150 # apply the modifiers
151 #pos = pos+pos_modifier
152 #readlen = readlen - readlen_modifier
153 if readlen in master_read_dict[tran]["ambig"]:
154 if pos in master_read_dict[tran]["ambig"][readlen]:
155 master_read_dict[tran]["ambig"][readlen][pos] += 1
156 else:
157 master_read_dict[tran]["ambig"][readlen][pos] = 1
158 else:
159 master_read_dict[tran]["ambig"][readlen] = {pos:1}
160 return ambiguously_mapped_reads
161
162
163 def get_mismatch_pos(md_tag,pos,readlen,master_read_dict,tran,readseq):
164 nucs = ["A","T","G","C"]
165 mismatches = {}
166 total_so_far = 0
167 prev_char = ""
168 for char in md_tag:
169 if char in nucs:
170 if prev_char != "":
171 total_so_far += int(prev_char)
172 prev_char = ""
173 mismatches[total_so_far+len(mismatches)] = (readseq[total_so_far+len(mismatches)])
174 else:
175 if char != "^" and char != "N":
176 if prev_char == "":
177 prev_char = char
178 else:
179 total_so_far += int(prev_char+char)
180 prev_char = ""
181 readlen_modifier = 0
182 pos_modifier = 0
183 five_ok = False
184 three_ok = False
185 while five_ok == False:
186 for i in range(0,readlen):
187 if i in mismatches:
188 pos_modifier += 1
189 readlen_modifier += 1
190 else:
191 five_ok = True
192 break
193 five_ok = True
194
195
196 while three_ok == False:
197 for i in range(readlen-1,0,-1):
198 if i in mismatches:
199 readlen_modifier += 1
200 else:
201 three_ok = True
202 break
203 three_ok = True
204
205
206 return (pos_modifier, readlen_modifier, mismatches)
207
208
209
210 def process_bam(bam_filepath, transcriptome_info_dict_path,outputfile):
211 desc = "NULL"
212 start_time = time.time()
213 study_dict ={}
214 nuc_count_dict = {"mapped":{},"unmapped":{}}
215 dinuc_count_dict = {}
216 threeprime_nuc_count_dict = {"mapped":{},"unmapped":{}}
217 read_length_dict = {}
218 unambig_read_length_dict = {}
219 unmapped_dict = {}
220 master_dict = {"unambiguous_non_coding_count":0,"unambiguous_coding_count":0,"current_dir":os.getcwd()}
221
222 transcriptome_info_dict = {}
223 connection = sqlite3.connect(transcriptome_info_dict_path)
224 cursor = connection.cursor()
225 cursor.execute("SELECT transcript,cds_start,cds_stop,length,strand,chrom,tran_type from transcripts;")
226 result = cursor.fetchall()
227 for row in result:
228 transcriptome_info_dict[str(row[0])] = {"cds_start":row[1],"cds_stop":row[2],"length":row[3],"strand":row[4],"chrom":row[5],"exons":[],"tran_type":row[6]}
229 #print list(transcriptome_info_dict)[:10]
230
231 cursor.execute("SELECT * from exons;")
232 result = cursor.fetchall()
233 for row in result:
234 transcriptome_info_dict[str(row[0])]["exons"].append((row[1],row[2]))
235
236 #it might be the case that there are no multimappers, so set this to 0 first to avoid an error, it will be overwritten later if there is multimappers
237 multimappers = 0
238 unmapped_reads = 0
239 unambiguous_coding_count = 0
240 unambiguous_non_coding_count = 0
241 trip_periodicity_reads = 0
242
243 final_offsets = {"fiveprime":{"offsets":{}, "read_scores":{}}, "threeprime":{"offsets":{}, "read_scores":{}}}
244 master_read_dict = {}
245 prev_seq = ""
246 process_chunk = {"read_name":[["placeholder_tran","1","28"]]}
247 mapped_reads = 0
248 ambiguously_mapped_reads = 0
249 master_trip_dict = {"fiveprime":{}, "threeprime":{}}
250 master_offset_dict = {"fiveprime":{}, "threeprime":{}}
251 master_metagene_stop_dict = {"fiveprime":{}, "threeprime":{}}
252
253
254 os.system(f'samtools sort -n {bam_filepath} -o {bam_filepath}_n_sorted.bam')
255
256 pysam.set_verbosity(0)
257 infile = pysam.Samfile(f"{bam_filepath}_n_sorted.bam", "rb")
258 header = infile.header["HD"]
259
260 unsorted = False
261 if "SO" in header:
262 print("Sorting order: "+header["SO"])
263 if header["SO"] != "queryname":
264 print("Sorting order is not queryname")
265 unsorted = True
266 else:
267 unsorted = True
268 if unsorted == True:
269 print ("ERROR: Bam file appears to be unsorted or not sorted by read name. To sort by read name use the command: samtools sort -n input.bam output.bam")
270 print (header,bam_filepath)
271 sys.exit()
272 total_bam_lines = 0
273 all_ref_ids = infile.references
274
275 for read in infile.fetch(until_eof=True):
276 total_bam_lines += 1
277 if not read.is_unmapped:
278 ref = read.reference_id
279 tran = (all_ref_ids[ref]).split(".")[0]
280 mapped_reads += 1
281 if mapped_reads%1000000 == 0:
282 print ("{} reads parsed at {}".format(mapped_reads,(time.time()-start_time)))
283 pos = read.reference_start
284 readname = read.query_name
285 read_tags = read.tags
286 if readname == list(process_chunk)[0]:
287 process_chunk[readname].append([tran,pos,read_tags])
288 #if the current read is different from previous reads send 'process_chunk' to the 'processor' function, then start 'process_chunk' over using current read
289 else:
290 if list(process_chunk)[0] != "read_name":
291
292 #At this point we work out readseq, we do this for multiple reasons, firstly so we don't count the sequence from a read multiple times, just because
293 # it aligns multiple times and secondly we only call read.seq once (read.seq is computationally expensive)
294 seq = read.seq
295 readlen = len(seq)
296
297 # Note if a read maps ambiguously it will still be counted toward the read length distribution (however it will only be counted once, not each time it maps)
298 if readlen not in read_length_dict:
299 read_length_dict[readlen] = 0
300 read_length_dict[readlen] += 1
301
302 if readlen not in nuc_count_dict["mapped"]:
303 nuc_count_dict["mapped"][readlen] = {}
304 if readlen not in threeprime_nuc_count_dict["mapped"]:
305 threeprime_nuc_count_dict["mapped"][readlen] = {}
306 if readlen not in dinuc_count_dict:
307 dinuc_count_dict[readlen] = {"AA":0, "TA":0, "GA":0, "CA":0,
308 "AT":0, "TT":0, "GT":0, "CT":0,
309 "AG":0, "TG":0, "GG":0, "CG":0,
310 "AC":0, "TC":0, "GC":0, "CC":0}
311
312 for i in range(0,len(seq)):
313 if i not in nuc_count_dict["mapped"][readlen]:
314 nuc_count_dict["mapped"][readlen][i] = {"A":0, "T":0, "G":0, "C":0, "N":0}
315 nuc_count_dict["mapped"][readlen][i][seq[i]] += 1
316
317 for i in range(0,len(seq)):
318 try:
319 dinuc_count_dict[readlen][seq[i:i+2]] += 1
320 except:
321 pass
322
323 for i in range(len(seq),0,-1):
324 dist = i-len(seq)
325 if dist not in threeprime_nuc_count_dict["mapped"][readlen]:
326 threeprime_nuc_count_dict["mapped"][readlen][dist] = {"A":0, "T":0, "G":0, "C":0, "N":0}
327 threeprime_nuc_count_dict["mapped"][readlen][dist][seq[dist]] += 1
328 ambiguously_mapped_reads += processor(process_chunk, master_read_dict, transcriptome_info_dict,master_dict,prev_seq, unambig_read_length_dict)
329 process_chunk = {readname:[[tran, pos, read_tags]]}
330 prev_seq = read.seq
331 else:
332 unmapped_reads += 1
333
334 # Add this unmapped read to unmapped_dict so we can see what the most frequent unmapped read is.
335 seq = read.seq
336 readlen = len(seq)
337 if seq in unmapped_dict:
338 unmapped_dict[seq] += 1
339 else:
340 unmapped_dict[seq] = 1
341
342 # Populate the nuc_count_dict with this unmapped read
343 if readlen not in nuc_count_dict["unmapped"]:
344 nuc_count_dict["unmapped"][readlen] = {}
345 for i in range(0,len(seq)):
346 if i not in nuc_count_dict["unmapped"][readlen]:
347 nuc_count_dict["unmapped"][readlen][i] = {"A":0, "T":0, "G":0, "C":0, "N":0}
348 nuc_count_dict["unmapped"][readlen][i][seq[i]] += 1
349
350 if readlen not in threeprime_nuc_count_dict["unmapped"]:
351 threeprime_nuc_count_dict["unmapped"][readlen] = {}
352
353 for i in range(len(seq),0,-1):
354 dist = i-len(seq)
355 if dist not in threeprime_nuc_count_dict["unmapped"][readlen]:
356 threeprime_nuc_count_dict["unmapped"][readlen][dist] = {"A":0, "T":0, "G":0, "C":0, "N":0}
357 threeprime_nuc_count_dict["unmapped"][readlen][dist][seq[dist]] += 1
358
359 #add stats about mapped/unmapped reads to file dict which will be used for the final report
360 master_dict["total_bam_lines"] = total_bam_lines
361 master_dict["mapped_reads"] = mapped_reads
362 master_dict["unmapped_reads"] = unmapped_reads
363 master_dict["ambiguously_mapped_reads"] = ambiguously_mapped_reads
364
365 if "read_name" in master_read_dict:
366 del master_read_dict["read_name"]
367 print ("BAM file processed")
368 print ("Creating metagenes, triplet periodicity plots, etc.")
369
370 for tran in master_read_dict:
371 try:
372 cds_start = int(0 if transcriptome_info_dict[tran]["cds_start"] is None else transcriptome_info_dict[tran]["cds_start"])
373 cds_stop = int(0 if transcriptome_info_dict[tran]["cds_stop"] is None else transcriptome_info_dict[tran]["cds_stop"])
374 # print(tran, type(cds_start))
375 except:
376 print("Exception: ", tran)
377 continue
378
379 tranlen = transcriptome_info_dict[tran]["length"]
380 #Use this to discard transcripts with no 5' leader or 3' trailer
381 if cds_start > 1 and cds_stop < tranlen and transcriptome_info_dict[tran]["tran_type"] == 1:
382 for primetype in ["fiveprime", "threeprime"]:
383 # Create the triplet periodicity and metainfo plots based on both the 5' and 3' ends of reads
384 for readlength in master_read_dict[tran]["unambig"]:
385 #print "readlength", readlength
386 # for each fiveprime postion for this readlength within this transcript
387 for raw_pos in master_read_dict[tran]["unambig"][readlength]:
388 #print "raw pos", raw_pos
389 trip_periodicity_reads += 1
390 if primetype == "fiveprime":
391 # get the five prime postion minus the cds start postion
392 real_pos = raw_pos-cds_start
393 rel_stop_pos = raw_pos-cds_stop
394 elif primetype == "threeprime":
395 real_pos = (raw_pos+readlength)-cds_start
396 rel_stop_pos = (raw_pos+readlength)-cds_stop
397 #get the readcount at the raw postion
398 readcount = master_read_dict[tran]["unambig"][readlength][raw_pos]
399 #print "readcount", readcount
400 frame = (real_pos%3)
401 if real_pos >= cds_start and real_pos <= cds_stop:
402 if readlength in master_trip_dict[primetype]:
403 master_trip_dict[primetype][readlength][str(frame)] += readcount
404 else:
405 master_trip_dict[primetype][readlength]= {"0":0.0,"1":0.0,"2":0.0}
406 master_trip_dict[primetype][readlength][str(frame)] += readcount
407 # now populate offset dict with the 'real_positions' upstream of cds_start, these will be used for metainfo dict
408 if real_pos > (-600) and real_pos < (601):
409 if readlength in master_offset_dict[primetype]:
410 if real_pos in master_offset_dict[primetype][readlength]:
411 #print "real pos in offset dict"
412 master_offset_dict[primetype][readlength][real_pos] += readcount
413 else:
414 #print "real pos not in offset dict"
415 master_offset_dict[primetype][readlength][real_pos] = readcount
416 else:
417 #initiliase with zero to avoid missing neighbours below
418 #print "initialising with zeros"
419 master_offset_dict[primetype][readlength]= {}
420 for i in range(-600,601):
421 master_offset_dict[primetype][readlength][i] = 0
422 master_offset_dict[primetype][readlength][real_pos] += readcount
423
424 # now populate offset dict with the 'real_positions' upstream of cds_start, these will be used for metainfo dict
425 if rel_stop_pos > (-600) and rel_stop_pos < (601):
426 if readlength in master_metagene_stop_dict[primetype]:
427 if rel_stop_pos in master_metagene_stop_dict[primetype][readlength]:
428 master_metagene_stop_dict[primetype][readlength][rel_stop_pos] += readcount
429 else:
430 master_metagene_stop_dict[primetype][readlength][rel_stop_pos] = readcount
431 else:
432 #initiliase with zero to avoid missing neighbours below
433 master_metagene_stop_dict[primetype][readlength] = {}
434 for i in range(-600,601):
435 master_metagene_stop_dict[primetype][readlength][i] = 0
436 master_metagene_stop_dict[primetype][readlength][rel_stop_pos] += readcount
437
438 # master trip dict is now made up of readlengths with 3 frames and a count associated with each frame
439 # create a 'score' for each readlength by putting the max frame count over the second highest frame count
440 for primetype in ["fiveprime", "threeprime"]:
441 for subreadlength in master_trip_dict[primetype]:
442 maxcount = 0
443 secondmaxcount = 0
444 for frame in master_trip_dict[primetype][subreadlength]:
445 if master_trip_dict[primetype][subreadlength][frame] > maxcount:
446 maxcount = master_trip_dict[primetype][subreadlength][frame]
447 for frame in master_trip_dict[primetype][subreadlength]:
448 if master_trip_dict[primetype][subreadlength][frame] > secondmaxcount and master_trip_dict[primetype][subreadlength][frame] != maxcount:
449 secondmaxcount = master_trip_dict[primetype][subreadlength][frame]
450 # a perfect score would be 0 meaning there is only a single peak, the worst score would be 1 meaning two highest peaks are the same height
451 master_trip_dict[primetype][subreadlength]["score"] = float(secondmaxcount)/float(maxcount)
452 #This part is to determine what offsets to give each read length
453 print ("Calculating offsets")
454 for primetype in ["fiveprime", "threeprime"]:
455 for readlen in master_offset_dict[primetype]:
456 accepted_len = False
457 max_relative_pos = 0
458 max_relative_count = 0
459 for relative_pos in master_offset_dict[primetype][readlen]:
460 # This line is to ensure we don't choose an offset greater than the readlength (in cases of a large peak far up/downstream)
461 if abs(relative_pos) < 10 or abs(relative_pos) > (readlen-10):
462 continue
463 if master_offset_dict[primetype][readlen][relative_pos] > max_relative_count:
464 max_relative_pos = relative_pos
465 max_relative_count = master_offset_dict[primetype][readlen][relative_pos]
466 #print "for readlen {} the max_relative pos is {}".format(readlen, max_relative_pos)
467 if primetype == "fiveprime":
468 # -3 to get from p-site to a-site, +1 to account for 1 based co-ordinates, resulting in -2 overall
469 final_offsets[primetype]["offsets"][readlen] = abs(max_relative_pos-2)
470 elif primetype == "threeprime":
471 # +3 to get from p-site to a-site, -1 to account for 1 based co-ordinates, resulting in +2 overall
472 final_offsets[primetype]["offsets"][readlen] = (max_relative_pos*(-1))+2
473 #If there are no reads in CDS regions for a specific length, it may not be present in master_trip_dict
474 if readlen in master_trip_dict[primetype]:
475 final_offsets[primetype]["read_scores"][readlen] = master_trip_dict[primetype][readlen]["score"]
476 else:
477 final_offsets[primetype]["read_scores"][readlen] = 0.0
478
479
480 master_read_dict["unmapped_reads"] = unmapped_reads
481 master_read_dict["offsets"] = final_offsets
482 master_read_dict["trip_periodicity"] = master_trip_dict
483 master_read_dict["desc"] = "Null"
484 master_read_dict["mapped_reads"] = mapped_reads
485 master_read_dict["nuc_counts"] = nuc_count_dict
486 master_read_dict["dinuc_counts"] = dinuc_count_dict
487 master_read_dict["threeprime_nuc_counts"] = threeprime_nuc_count_dict
488 master_read_dict["metagene_counts"] = master_offset_dict
489 master_read_dict["stop_metagene_counts"] = master_metagene_stop_dict
490 master_read_dict["read_lengths"] = read_length_dict
491 master_read_dict["unambig_read_lengths"] = unambig_read_length_dict
492 master_read_dict["coding_counts"] = master_dict["unambiguous_coding_count"]
493 master_read_dict["noncoding_counts"] = master_dict["unambiguous_non_coding_count"]
494 master_read_dict["ambiguous_counts"] = master_dict["ambiguously_mapped_reads"]
495 master_read_dict["frequent_unmapped_reads"] = (sorted(unmapped_dict.items(), key=operator.itemgetter(1)))[-2000:]
496 master_read_dict["cutadapt_removed"] = 0
497 master_read_dict["rrna_removed"] = 0
498 #If no reads are removed by minus m there won't be an entry in the log file, so initiliase with 0 first and change if there is a line
499 master_read_dict["removed_minus_m"] = 0
500 master_dict["removed_minus_m"] = 0
501 # We work out the total counts for 5', cds 3' for differential translation here, would be better to do thisn in processor but need the offsets
502 master_read_dict["unambiguous_all_totals"] = {}
503 master_read_dict["unambiguous_fiveprime_totals"] = {}
504 master_read_dict["unambiguous_cds_totals"] = {}
505 master_read_dict["unambiguous_threeprime_totals"] = {}
506
507 master_read_dict["ambiguous_all_totals"] = {}
508 master_read_dict["ambiguous_fiveprime_totals"] = {}
509 master_read_dict["ambiguous_cds_totals"] = {}
510 master_read_dict["ambiguous_threeprime_totals"] = {}
511 print ("calculating transcript counts")
512 for tran in master_read_dict:
513 if tran in transcriptome_info_dict:
514 five_total = 0
515 cds_total = 0
516 three_total = 0
517
518 ambig_five_total = 0
519 ambig_cds_total = 0
520 ambig_three_total = 0
521
522 cds_start = transcriptome_info_dict[tran]["cds_start"]
523 cds_stop = transcriptome_info_dict[tran]["cds_stop"]
524
525 for readlen in master_read_dict[tran]["unambig"]:
526 if readlen in final_offsets["fiveprime"]["offsets"]:
527 offset = final_offsets["fiveprime"]["offsets"][readlen]
528 else:
529 offset = 15
530 for pos in master_read_dict[tran]["unambig"][readlen]:
531 real_pos = pos+offset
532 if cds_start is None or cds_stop is None:
533 three_total += master_read_dict[tran]["unambig"][readlen][pos]
534 else:
535 if real_pos <cds_start:
536 five_total += master_read_dict[tran]["unambig"][readlen][pos]
537 elif real_pos >=cds_start and real_pos <= cds_stop:
538 cds_total += master_read_dict[tran]["unambig"][readlen][pos]
539 elif real_pos > cds_stop:
540 three_total += master_read_dict[tran]["unambig"][readlen][pos]
541 master_read_dict["unambiguous_all_totals"][tran] = five_total+cds_total+three_total
542 master_read_dict["unambiguous_fiveprime_totals"][tran] = five_total
543 master_read_dict["unambiguous_cds_totals"][tran] = cds_total
544 master_read_dict["unambiguous_threeprime_totals"][tran] = three_total
545
546 for readlen in master_read_dict[tran]["ambig"]:
547 if readlen in final_offsets["fiveprime"]["offsets"]:
548 offset = final_offsets["fiveprime"]["offsets"][readlen]
549 else:
550 offset = 15
551 for pos in master_read_dict[tran]["ambig"][readlen]:
552 if cds_start is None or cds_stop is None:
553 ambig_three_total += master_read_dict[tran]["ambig"][readlen][pos]
554 else:
555 real_pos = pos+offset
556 if real_pos < cds_start:
557 ambig_five_total += master_read_dict[tran]["ambig"][readlen][pos]
558 elif real_pos >=cds_start and real_pos <= cds_stop:
559 ambig_cds_total += master_read_dict[tran]["ambig"][readlen][pos]
560 elif real_pos > cds_stop:
561 ambig_three_total += master_read_dict[tran]["ambig"][readlen][pos]
562
563 master_read_dict["ambiguous_all_totals"][tran] = five_total+cds_total+three_total+ambig_five_total+ambig_cds_total+ambig_three_total
564 master_read_dict["ambiguous_fiveprime_totals"][tran] = five_total+ambig_five_total
565 master_read_dict["ambiguous_cds_totals"][tran] = cds_total+ambig_cds_total
566 master_read_dict["ambiguous_threeprime_totals"][tran] = three_total+ambig_three_total
567
568 print ("Writing out to sqlite file")
569 sqlite_db = SqliteDict(outputfile,autocommit=False)
570 for key in master_read_dict:
571 sqlite_db[key] = master_read_dict[key]
572 sqlite_db["description"] = desc
573 sqlite_db.commit()
574 sqlite_db.close()
575
576
577 if __name__ == "__main__":
578 if len(sys.argv) <= 2:
579 print ("Usage: python bam_to_sqlite.py <path_to_bam_file> <path_to_organism.sqlite> <file_description (optional)>")
580 sys.exit()
581 bam_filepath = sys.argv[1]
582 annotation_sqlite_filepath = sys.argv[2]
583 desc = sys.argv[3]
584 outputfile = sys.argv[4]
585 process_bam(bam_filepath,annotation_sqlite_filepath,outputfile)