Mercurial > repos > glogobyte > isoread

comparison mirbase_ultra_v2.py @ 0:99e8a03f8802 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author glogobyte
date Fri, 16 Oct 2020 18:53:05 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:99e8a03f8802
1 from mirbase_functions import *
2 from mirbase_graphs import *
3 import itertools
4 import time
5 import sys
6 import os
7 import urllib.request
8 import gzip
9 from multiprocessing import Process, Queue, Lock, Pool, Manager, Value
10 import subprocess
11 import argparse
12 from collections import OrderedDict
13 from matplotlib.backends.backend_pdf import PdfPages
14 import pandas as pd
15 from math import pi
16 import numpy as np
17 import matplotlib.pyplot as plt
18 from matplotlib.ticker import PercentFormatter
19 import seaborn as sns
20 import scipy.stats as stats
21 from plotnine import *
22 import math
23 import re
24 import matplotlib.ticker as mtick
25 import copy
26
27 subprocess.call(['mkdir','-p', 'split1','split2','split3','split4','split11','split12','Counts','Diff/temp_con','Diff/temp_tre','Diff/n_temp_con','Diff/n_temp_tre'])
28
29 parser = argparse.ArgumentParser()
30 parser.add_argument("-analysis", "--anal", help="choose type of analysis", action="store")
31 parser.add_argument("-con", "--control", help="input fastq file", nargs='+', default=[])
32 parser.add_argument("-tre", "--treated", help="input fastq file", nargs='+', default=[] )
33 parser.add_argument("-tool_dir", "--tool_directory", help="tool directory path", action="store")
34 parser.add_argument("-gen", "--org_name", help="tool directory path", action="store")
35 parser.add_argument("-program", "--pro", help="choose type of analysis", action="store")
36 parser.add_argument("-f", "--flag", help="choose the database", action="store")
37 parser.add_argument("-umis", "--umi", help="choose the database", action="store")
38 parser.add_argument("-percentage", "--per", help="choose the database", action="store")
39 parser.add_argument("-counts", "--count", help="choose the database", action="store")
40 parser.add_argument("-name1", "--n1", help="choose the database", action="store")
41 parser.add_argument("-name2", "--n2", help="choose the database", action="store")
42 args = parser.parse_args()
43
44 #########################################################################################################################################
45
46 def collapse_sam(path):
47
48 ini_sam=read(path,0)
49 main_sam = [x.rstrip("\n").split("\t") for x in ini_sam if "@" not in x.split("\t")[0]]
50 intro_sam = [x.rstrip("\n").split("\t") for x in ini_sam if "@" in x.split("\t")[0]]
51
52 uni_seq = []
53 for x in main_sam:
54
55 if [x[2], x[9]] not in uni_seq:
56 uni_seq.append([x[2], x[9]])
57
58 new_main_sam=[]
59 incr_num=0
60 for i in range(len(uni_seq)):
61 count=0
62 incr_num+=1
63 for y in main_sam:
64 if uni_seq[i][1]==y[9] and uni_seq[i][0]==y[2]:
65 count+=1
66 temp=y
67 temp[10]="~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
68 temp[0]=str(incr_num)+"-"+str(count)
69 new_main_sam.append(temp)
70
71 new_sam=intro_sam+new_main_sam
72
73 return new_sam
74
75 #################################################################################################################################################################################################################
76
77 def duplicate_chroms_isoforms(List):
78
79 dupes=[]
80
81 for num in range(len(List)):
82
83 if [List[num][9],List[num][0],List[num][2]] not in dupes :
84 dupes.append([List[num][9],List[num][0],List[num][2]])
85
86 for x in List:
87 for y in dupes:
88 if x[9]==y[0] and x[0]==y[1] and x[2].split("_")[0]==y[2].split("_")[0] and x[2]!=y[2]:
89 y.append(x[2])
90
91
92 double_List = [x[:] for x in List]
93
94 chr_order=[]
95 for x in dupes:
96 temp = []
97 for i in range(2,len(x)):
98 if x[i].split("chr")[1].split("(")[0].isdigit():
99 temp.append(int(x[i].split("chr")[1].split("(")[1][0]+x[i].split("chr")[1].split("(")[0]))
100 else:
101 temp.append(x[i].split("chr")[1][0:4])
102
103 for z in temp:
104 if 'X(-)'==z or 'Y(-)'==z or 'X(+)'==z or 'Y(+)'==z:
105 temp = [str(j) for j in temp]
106 temp=list(set(temp))
107 temp.sort()
108 chr_order.append(temp)
109
110 final_dupes=[]
111 for i in range(len(dupes)):
112 final_dupes.append([dupes[i][0],dupes[i][2].split("_")[0],dupes[i][1]])
113 for x in chr_order[i]:
114 result = re.match("[-+]?\d+$", str(x))
115 if len(chr_order[i]) == len(set(chr_order[i])):
116 if result is not None:
117
118 if int(x)<0:
119 final_dupes[i][1]=final_dupes[i][1]+"_chr"+str(abs(int(x)))+"(-)"
120 else:
121 final_dupes[i][1] = final_dupes[i][1] + "_chr" + str(abs(int(x)))+"(+)"
122 else:
123 final_dupes[i][1] = final_dupes[i][1] + "_chr" + str(x)
124 else:
125 if result is not None:
126 if int(x) < 0:
127 final_dupes[i][1] = final_dupes[i][1] +dupes[i][2].split("_")[1]+ "_chr" + str(abs(int(x))) + "(-)"
128 else:
129 final_dupes[i][1] = final_dupes[i][1] +dupes[i][2].split("_")[1]+ "_chr" + str(abs(int(x))) + "(+)"
130 else:
131 final_dupes[i][1] = final_dupes[i][1] +dupes[i][2].split("_")[1]+ "_chr" + str(x)
132
133 final_dupes.sort()
134 final_dupes=list(final_dupes for final_dupes,_ in itertools.groupby(final_dupes))
135
136 for i in range(len(double_List)):
137 for x in final_dupes:
138
139 if double_List[i][9] == x[0] and double_List[i][0] == x[2] and len(double_List[i][2].split("_")) >3 and double_List[i][2].split("_")[0]==x[1].split("_")[0]:
140 gg=str("_"+double_List[i][2].split("_")[-2]+"_"+double_List[i][2].split("_")[-1])
141 double_List[i][2] = x[1]+gg
142
143 if double_List[i][9]==x[0] and double_List[i][0]== x[2] and len(double_List[i][2].split("_"))==3 and double_List[i][2].split("_")[0]==x[1].split("_")[0]:
144 double_List[i][2]=x[1]
145 List[i][2] = x[1]
146
147 List.sort()
148 new_list=list(List for List,_ in itertools.groupby(List))
149
150 double_List.sort()
151 new_double_List = list(double_List for double_List, _ in itertools.groupby(double_List))
152
153 return new_list, new_double_List
154
155
156 #############################################################################################################################################################################################################
157
158 def sam(mature_mirnas,path,name,con,l,samples,data,names,unmap_seq,samples_mirna_names,deseq,LHE_names,umi,ini_sample,unmap_counts):
159
160 # read the sam file
161 ini_sam=read(path,0)
162 new_main_sam = [x.rstrip("\n").split("\t") for x in ini_sam if "@" not in x.split("\t")[0]]
163 unique_seq = [x for x in new_main_sam if x[1] == '0' and len(x[9])>=18 and len(x[9])<=26]
164
165 sorted_uni_arms = []
166
167 for i in range(len(mature_mirnas)):
168 tmp_count_reads = 0 # calculate the total number of reads
169 tmp_count_seq = 0 # calculate the total number of sequences
170 for j in range(len(unique_seq)):
171
172 if "{" in unique_seq[j][2].split("_")[0]:
173 official=unique_seq[j][2].split("_")[0][:-4]
174 else:
175 official=unique_seq[j][2].split("_")[0]
176
177 if mature_mirnas[i].split(" ")[0][1:] == official:
178
179 temp_mature = mature_mirnas[i+1].strip().replace("U", "T")
180 off_part = longestSubstring(temp_mature, unique_seq[j][9])
181
182 mat_diff = temp_mature.split(off_part)
183 mat_diff = [len(mat_diff[0]), len(mat_diff[1])]
184
185 unique_diff = unique_seq[j][9].split(off_part)
186 unique_diff = [len(unique_diff[0]), len(unique_diff[1])]
187
188 # Problem with hsa-miR-8485
189 if mat_diff[1]!=0 and unique_diff[1]!=0:
190 unique_seq[j]=1
191 pre_pos = 0
192 post_pos = 0
193
194 elif mat_diff[0]!=0 and unique_diff[0]!=0:
195 unique_seq[j]=1
196 pre_pos = 0
197 post_pos = 0
198
199 else:
200 pre_pos = mat_diff[0]-unique_diff[0]
201 post_pos = unique_diff[1]-mat_diff[1]
202 tmp_count_reads = tmp_count_reads + int(unique_seq[j][0].split("-")[1])
203 tmp_count_seq = tmp_count_seq+1
204
205 if pre_pos != 0 or post_pos != 0:
206 if pre_pos == 0:
207 unique_seq[j][2] = unique_seq[j][2] + "_" +str(pre_pos) + "_" + '{:+d}'.format(post_pos)
208 elif post_pos == 0:
209 unique_seq[j][2] = unique_seq[j][2] + "_" + '{:+d}'.format(pre_pos) + "_" + str(post_pos)
210 else:
211 unique_seq[j][2] = unique_seq[j][2]+"_"+'{:+d}'.format(pre_pos)+"_"+'{:+d}'.format(post_pos)
212
213 for x in range(unique_seq.count(1)):
214 unique_seq.remove(1)
215 if tmp_count_reads != 0 and tmp_count_seq != 0:
216 sorted_uni_arms.append([mature_mirnas[i].split(" ")[0][1:], tmp_count_seq, tmp_count_reads])
217 sorted_uni_arms = sorted(sorted_uni_arms, key=lambda x: x[1], reverse=True)
218 dedup_unique_seq,double_fil_uni_seq=duplicate_chroms_isoforms(unique_seq)
219
220 for y in sorted_uni_arms:
221 counts=0
222 seqs=0
223 for x in double_fil_uni_seq:
224 if y[0]==x[2].split("_")[0]:
225 counts+=int(x[0].split("-")[1])
226 seqs+=1
227
228 y[1]=seqs
229 y[2]=counts
230
231 LHE=[]
232 l.acquire()
233 if con=="c":
234 LHE.extend(z[2] for z in double_fil_uni_seq)
235 for y in double_fil_uni_seq:
236 samples_mirna_names.append([y[2],y[9]])
237 deseq.append([[x[2], x[0].split('-')[1], x[9]] for x in double_fil_uni_seq])
238 LHE_names.extend(LHE)
239 unmap_seq.value += sum([1 for x in new_main_sam if x[1] == '4'])
240 unmap_counts.value += sum([int(x[0].split("-")[1]) for x in new_main_sam if x[1] == '4'])
241 names.append(name)
242 samples.append(dedup_unique_seq)
243 data.append([con,name,double_fil_uni_seq,sorted_uni_arms])
244 ini_sample.append(new_main_sam)
245
246 if con=="t":
247 LHE.extend(z[2] for z in double_fil_uni_seq)
248 for y in double_fil_uni_seq:
249 samples_mirna_names.append([y[2],y[9]])
250 deseq.append([[x[2], x[0].split('-')[1], x[9]] for x in double_fil_uni_seq])
251 LHE_names.extend(LHE)
252 unmap_seq.value += sum([1 for x in new_main_sam if x[1] == '4'])
253 unmap_counts.value += sum([int(x[0].split("-")[1]) for x in new_main_sam if x[1] == '4'])
254 names.append(name)
255 samples.append(dedup_unique_seq)
256 data.append([con,name,double_fil_uni_seq,sorted_uni_arms])
257 ini_sample.append(new_main_sam)
258 l.release()
259
260
261 ######################################################################################################################################
262 """
263
264 Read a sam file from Bowtie and do the followings:
265
266 1) Remove reverse stranded mapped reads
267 2) Remove unmapped reads
268 3) Remove all sequences with reads less than 11 reads
269 4) Sort the arms with the most sequences in decreading rate
270 5) Sort the sequences of every arm with the most reads in decreasing rate
271 6) Calculate total number of sequences of every arm
272 7) Calculate total number of reads of sequences of every arm.
273 8) Store all the informations in a txt file
274
275 """
276
277 def non_sam(mature_mirnas,path,name,con,l,data,names,n_deseq,n_samples_mirna_names,n_LHE_names):
278
279 ini_sam=read(path,0)
280 new_main_sam = [x.rstrip("\n").split("\t") for x in ini_sam if "@" not in x.split("\t")[0]]
281 unique_seq=[]
282 unique_seq = [x for x in new_main_sam if x[1] == '4' and len(x[9])>=18 and len(x[9])<=26]
283
284 uni_seq=[]
285 # Calculate the shifted positions for every isomir and add them to the name of it
286 sorted_uni_arms = []
287 for i in range(1,len(mature_mirnas),2):
288 tmp_count_reads = 0 # calculate the total number of reads
289 tmp_count_seq = 0 # calculate the total number of sequences
290
291 for j in range(len(unique_seq)):
292
293 temp_mature = mature_mirnas[i].strip().replace("U", "T")
294
295 if temp_mature in unique_seq[j][9]:
296
297 off_part = longestSubstring(temp_mature, unique_seq[j][9])
298
299 mat_diff = temp_mature.split(off_part)
300 mat_diff = [len(mat_diff[0]), len(mat_diff[1])]
301
302 unique_diff = unique_seq[j][9].split(off_part)
303 if len(unique_diff)<=2:
304 unique_diff = [len(unique_diff[0]), len(unique_diff[1])]
305
306 pre_pos = mat_diff[0]-unique_diff[0]
307 post_pos = unique_diff[1]-mat_diff[1]
308
309 lengthofmir = len(off_part) + post_pos
310 if pre_pos == 0:
311 tmp_count_reads = tmp_count_reads + int(unique_seq[j][0].split("-")[1])
312 tmp_count_seq = tmp_count_seq + 1
313
314 if pre_pos == 0:
315
316 t_name=unique_seq[j].copy()
317 t_name[2]=mature_mirnas[i - 1].split(" ")[0][1:] + "__" + str(pre_pos) + "_" + '{:+d}'.format(post_pos) + "_" + str(unique_seq[j][9][len(off_part):])
318 uni_seq.append(t_name)
319
320
321 if tmp_count_reads != 0 and tmp_count_seq != 0:
322 sorted_uni_arms.append([mature_mirnas[i-1].split(" ")[0][1:], tmp_count_seq, tmp_count_reads])
323
324
325 sorted_uni_arms = sorted(sorted_uni_arms, key=lambda x: x[1], reverse=True)
326 unique_seq = list(map(list, OrderedDict.fromkeys(map(tuple,uni_seq))))
327
328 LHE=[]
329
330 l.acquire()
331 if con=="c":
332 LHE.extend(x[2] for x in unique_seq if x[2]!="*")
333 for x in unique_seq:
334 if x[2]!="*":
335 n_samples_mirna_names.append([x[2],x[9]])
336 n_deseq.append([[x[2], x[0].split('-')[1], x[9]] for x in unique_seq if x[2]!="*"])
337 n_LHE_names.extend(LHE)
338 names.append(name)
339 data.append([con,name,unique_seq,sorted_uni_arms])
340
341
342 if con=="t":
343 LHE.extend(x[2] for x in unique_seq if x[2]!="*")
344 for x in unique_seq:
345 if x[2]!="*":
346 n_samples_mirna_names.append([x[2],x[9]])
347 n_deseq.append([[x[2], x[0].split('-')[1], x[9]] for x in unique_seq if x[2]!="*"])
348 n_LHE_names.extend(LHE)
349 names.append(name)
350 data.append([con,name,unique_seq,sorted_uni_arms])
351 l.release()
352
353 #####################################################################################################################################################################################################################
354 def deseq2_temp(samples_mirna_names,deseq,con,l):
355
356 samples_mirna_names.sort(key=lambda x:[0])
357 for i in range(len(deseq)):
358 for y in samples_mirna_names:
359 flag = 0
360 for x in deseq[i]:
361 if y[0] == x[0]:
362 flag = 1
363 break
364
365 if flag == 0:
366 deseq[i].append([y[0], "0", y[1]])
367
368 [deseq[i].sort(key=lambda x: x[0]) for i, _ in enumerate(deseq)]
369 deseq_final = [[x[0],x[2]] for x in deseq[0]]
370 [deseq_final[z].append(deseq[i][j][1]) for z,_ in enumerate(deseq_final) for i, _ in enumerate(deseq) for j,_ in enumerate(deseq[i]) if deseq_final[z][0] == deseq[i][j][0]]
371
372 l.acquire()
373 if con=="c":
374 q1.put(deseq_final)
375
376 if con=="t":
377 q2.put(deseq_final)
378 l.release()
379
380
381 ####################################################################################################################################################################################################################
382
383 def main_temp(LH2E, LH2E_names, LH8E, LH8E_names,flag,names_con,names_tre,filter_LH8E,filter_LH2E,raw_LH8E,raw_LH2E):
384
385 LH8E_add_names = [x for x in LH2E_names if x not in LH8E_names]
386 LH2E_add_names = [x for x in LH8E_names if x not in LH2E_names]
387
388 LH8E_add_names.sort()
389 LH2E_add_names.sort()
390 LH8E_add_names = list(LH8E_add_names for LH8E_add_names,_ in itertools.groupby(LH8E_add_names))
391 LH2E_add_names = list(LH2E_add_names for LH2E_add_names,_ in itertools.groupby(LH2E_add_names))
392
393 LH2E.sort()
394 LH8E.sort()
395 LH2E = list(LH2E for LH2E,_ in itertools.groupby(LH2E))
396 LH8E = list(LH8E for LH8E,_ in itertools.groupby(LH8E))
397
398 print("LHE_names")
399 print([len(LH8E_add_names),len(LH2E_add_names)])
400 print([len(LH8E),len(LH2E)])
401
402 zeros=["0"]*(len(LH8E[0])-2)
403 [LH8E_add_names[i].extend(zeros) for i,_ in enumerate(LH8E_add_names)]
404 LH8E=LH8E+LH8E_add_names
405
406 zeros=["0"]*(len(LH2E[0])-2)
407 [LH2E_add_names[i].extend(zeros) for i,_ in enumerate(LH2E_add_names)]
408 LH2E=LH2E+LH2E_add_names
409
410 dupes=[]
411 final_LH2E =[]
412
413 for num,_ in enumerate(LH2E):
414
415 if LH2E[num][1] not in final_LH2E and LH2E[num][0] not in final_LH2E:
416 final_LH2E.append(LH2E[num][1])
417 final_LH2E.append(LH2E[num][0])
418 else:
419 dupes.append(LH2E[num][1])
420
421
422 dupes=list(set(dupes))
423
424 dupes=[[x] for x in dupes]
425
426 for x in LH2E:
427 for y in dupes:
428 if x[1]==y[0]:
429 fl=0
430 if len(y)==1:
431 y.append(x[0])
432 else:
433 for i in range(1,len(y)):
434 if y[i].split("_")[0]==x[0].split("_")[0]:
435 fl=1
436 if len(x[0])<len(y[i]):
437 del y[i]
438 y.append(x[0])
439 break
440
441 if fl==0:
442 y.append((x[0]))
443
444 for y in dupes:
445 if len(y)>2:
446 for i in range(len(y)-1,1,-1):
447 y[1]=y[1]+"/"+y[i]
448 del y[i]
449
450 for x in LH2E:
451 for y in dupes:
452 if x[1]==y[0]:
453 x[0]=y[1]
454
455 for x in LH8E:
456 for y in dupes:
457 if x[1]==y[0]:
458 x[0]=y[1]
459
460
461 LH2E.sort()
462 LH2E=list(LH2E for LH2E,_ in itertools.groupby(LH2E))
463
464 LH8E.sort()
465 LH8E=list(LH8E for LH8E,_ in itertools.groupby(LH8E))
466
467 LH8E_new=[]
468 LH2E_new=[]
469
470 if int(args.per)!=-1:
471 percent=int(args.per)/100
472 print(percent)
473 print(args.count)
474
475 c_col_filter=round(percent*(len(LH2E[1])-2))
476 t_col_filter=round(percent*(len(LH8E[1])-2))
477
478 for i, _ in enumerate(LH2E):
479 c_cols=0
480 t_cols=0
481
482 c_cols=sum([1 for j in range(len(LH2E[i])-2) if int(LH2E[i][j+2])>=int(args.count)])
483 t_cols=sum([1 for j in range(len(LH8E[i])-2) if int(LH8E[i][j+2])>=int(args.count)])
484
485 if c_cols>=c_col_filter or t_cols>=t_col_filter:
486 LH8E_new.append(LH8E[i])
487 LH2E_new.append(LH2E[i])
488
489 filter_LH2E.extend(LH2E_new)
490 filter_LH8E.extend(LH8E_new)
491 raw_LH2E.extend(LH2E)
492 raw_LH8E.extend(LH8E)
493
494 ##################################################################################################################################################################################################################
495
496 def write_main(raw_LH2E, raw_LH8E, fil_LH2E, fil_LH8E, names_con, names_tre, flag):
497
498 if flag == 1 and int(args.per)!=-1:
499 fp = open('Counts/Filtered '+args.n2 +' Templated Counts', 'w')
500 fp.write("Name\t")
501 fp.write("Sequence")
502 for y in names_tre:
503 fp.write("\t"+y)
504
505 for x in fil_LH8E:
506 fp.write("\n%s" % "\t".join(x))
507 fp.close()
508
509 fp = open('Counts/Filtered '+args.n1+' Templated Counts', 'w')
510 fp.write("Name\t")
511 fp.write("Sequence")
512 for y in names_con:
513 fp.write("\t"+y)
514
515 for x in fil_LH2E:
516 fp.write("\n%s" % "\t".join(x))
517 fp.close()
518
519
520 if flag == 2 and int(args.per)!=-1:
521 fp = open('Counts/Filtered '+args.n2+' Non-Templated Counts', 'w')
522 fp.write("Name\t")
523 fp.write("Sequence")
524 for y in names_tre:
525 fp.write("\t"+y)
526
527
528 for x in fil_LH8E:
529 fp.write("\n%s" % "\t".join(x))
530 fp.close()
531
532 fp = open('Counts/Filtered '+args.n1+' Non-Templated Counts', 'w')
533 fp.write("Name\t")
534 fp.write("Sequence")
535 for y in names_con:
536 fp.write("\t"+y)
537
538 for x in fil_LH2E:
539 fp.write("\n%s" % "\t".join(x))
540 fp.close()
541
542
543 if flag == 1:
544 fp = open('Counts/Raw '+args.n2+' Templated Counts', 'w')
545 fp.write("Name\t")
546 fp.write("Sequence")
547 for y in names_tre:
548 fp.write("\t"+y)
549
550 for x in raw_LH8E:
551 fp.write("\n%s" % "\t".join(x))
552 fp.close()
553
554 fp = open('Counts/Raw '+args.n1+' Templated Counts', 'w')
555 fp.write("Name\t")
556 fp.write("Sequence")
557 for y in names_con:
558 fp.write("\t"+y)
559
560 for x in raw_LH2E:
561 fp.write("\n%s" % "\t".join(x))
562 fp.close()
563
564
565 if flag == 2:
566 fp = open('Counts/Raw '+args.n2+' Non-Templated Counts', 'w')
567 fp.write("Name\t")
568 fp.write("Sequence")
569 for y in names_tre:
570 fp.write("\t"+y)
571
572
573 for x in raw_LH8E:
574 fp.write("\n%s" % "\t".join(x))
575 fp.close()
576
577 fp = open('Counts/Raw '+args.n1+' Non-Templated Counts', 'w')
578 fp.write("Name\t")
579 fp.write("Sequence")
580 for y in names_con:
581 fp.write("\t"+y)
582
583 for x in raw_LH2E:
584 fp.write("\n%s" % "\t".join(x))
585 fp.close()
586
587
588 #########################################################################################################################################
589
590 def ssamples(names,samp,folder,pro):
591
592 for i in range(2,len(samp[0])):
593
594 fp = open(folder+names[i-2]+'.txt','w')
595 fp.write("miRNA id"+"\t"+names[i-2]+"\n")
596
597 for x in samp:
598 fp.write("%s" % "\t".join([x[0],x[i]])+"\n")
599 fp.close()
600
601 ##################################################################################################################
602
603 def DB_write(con,name,unique_seq,sorted_uni_arms,f):
604
605 if f==1:
606 # Write a txt file with all the information
607 if con=="c":
608 fp = open('split1/'+name, 'w')
609
610 fp.write("%s\t%-42s\t%s\n\n" % ("Number of Reads","Name of isomir","Sequence"))
611 if con=="t":
612 fp = open('split2/'+name, 'w')
613 fp.write("%s\t%-42s\t%s\n\n" % ("Number of Reads","Name of isomir","Sequence"))
614
615
616 for i in range(len(sorted_uni_arms)):
617 temp = []
618 for j in range(len(unique_seq)):
619
620 if sorted_uni_arms[i][0] in unique_seq[j][2].split("_")[0]:
621
622 temp.append(unique_seq[j])
623
624 temp = sorted(temp, key=lambda x: int(x[0].split('-')[1]), reverse=True)
625 fp.write("*********************************************************************************************************\n")
626 fp.write("%-8s\t%-22s\t%-25s\t%-30s\t%s\n" % ("|",str(sorted_uni_arms[i][0]),"Sequence count = "+str(sorted_uni_arms[i][1]),"Total reads = "+str(sorted_uni_arms[i][2]),"|"))
627 fp.write("*********************************************************************************************************\n\n")
628 [fp.write("%-8s\t%-40s\t%s\n" % (x[0].split("-")[1], x[2],x[9])) for x in temp]
629 fp.write("\n" + "\n")
630 fp.close()
631
632 if f==2:
633
634 if con=="c":
635 fp = open('split3/'+name, 'w')
636 fp.write("%s\t%-42s\t%s\n\n" % ("Number of Reads","Name of isomir","Sequence"))
637 if con=="t":
638 fp = open('split4/'+name, 'w')
639 fp.write("%s\t%-42s\t%s\n\n" % ("Number of Reads","Name of isomir","Sequence"))
640
641
642 for i in range(len(sorted_uni_arms)):
643 temp = []
644 for j in range(len(unique_seq)):
645 if sorted_uni_arms[i][0]==unique_seq[j][2].split("__")[0]:
646 temp.append(unique_seq[j])
647 if temp!=[]:
648 temp = sorted(temp, key=lambda x: int(x[0].split('-')[1]), reverse=True)
649 fp.write("*********************************************************************************************************\n")
650 fp.write("%-8s\t%-22s\t%-25s\t%-30s\t%s\n" % ("|",str(sorted_uni_arms[i][0]),"Sequence count = "+str(sorted_uni_arms[i][1]),"Total reads = "+str(sorted_uni_arms[i][2]),"|"))
651 fp.write("*********************************************************************************************************\n\n")
652 [fp.write("%-8s\t%-40s\t%s\n" % (x[0].split("-")[1], x[2],x[9])) for x in temp]
653 fp.write("\n" + "\n")
654 fp.close()
655
656
657 ##########################################################################################################################
658
659 def new_mat_seq(pre_unique_seq,mat_mirnas,l):
660
661 unique_iso = []
662 for x in pre_unique_seq:
663 if len(x[2].split("_"))==3:
664 for y in pre_unique_seq:
665 if x[2] in y[2] and int(x[0].split("-")[1])<int(y[0].split("-")[1]):
666 if any(y[2] in lst2 for lst2 in unique_iso)==False:
667 y[2]=">"+y[2]
668 unique_iso.append(y)
669 l.acquire()
670 for x in unique_iso:
671 mat_mirnas.append(x[2])
672 mat_mirnas.append(x[9])
673 l.release()
674
675 #########################################################################################################################
676 def pie_non_temp(merge_LH2E,merge_non_LH2E,merge_LH8E,merge_non_LH8E,c_unmap,t_unmap,c_unmap_counts,t_unmap_counts):
677
678 c_samples=[[x[0],x[1],sum(int(i) for i in x[2:])] for x in merge_LH2E]
679 t_samples=[[x[0],x[1],sum(int(i) for i in x[2:])] for x in merge_LH8E]
680 c_non_samples=[[x[0],x[1],sum(int(i) for i in x[2:])] for x in merge_non_LH2E]
681 t_non_samples=[[x[0],x[1],sum(int(i) for i in x[2:])] for x in merge_non_LH8E]
682
683 c_templ = 0
684 c_tem_counts = 0
685 c_mature = 0
686 c_mat_counts = 0
687 t_templ = 0
688 t_tem_counts = 0
689 t_mature = 0
690 t_mat_counts = 0
691
692 c_non = len(c_non_samples)
693 c_non_counts = sum(x[2] for x in c_non_samples)
694 t_non = len(t_non_samples)
695 t_non_counts = sum(x[2] for x in t_non_samples)
696
697 c_unmap = c_unmap - c_non
698 t_unmap = c_unmap - t_non
699
700 c_unmap_counts=c_unmap_counts - c_non_counts
701 t_unmap_counts=t_unmap_counts - t_non_counts
702
703
704 for x in c_samples:
705
706 if "/" not in x[0]:
707 if "chr" in x[0].split("_")[-1]:
708 c_mature+=1
709 c_mat_counts += x[2]
710 else:
711 c_templ+=1
712 c_tem_counts += x[2]
713 else:
714 f=0
715 for y in x[0].split("/"):
716 if "chr" in y.split("_")[-1]:
717 c_mature+=1
718 c_mat_counts += x[2]
719 f=1
720 break
721 if f==0:
722 c_templ+=1
723 c_tem_counts += x[2]
724
725 for x in t_samples:
726
727 if "/" not in x[0]:
728 if "chr" in x[0].split("_")[-1]:
729 t_mature+=1
730 t_mat_counts += x[2]
731 else:
732 t_templ+=1
733 t_tem_counts += x[2]
734 else:
735 f=0
736 for y in x[0].split("/"):
737 if "chr" in y.split("_")[-1]:
738 t_mature+=1
739 t_mat_counts += x[2]
740 f=1
741 break
742 if f==0:
743 t_templ+=1
744 t_tem_counts += x[2]
745
746 fig = plt.figure(figsize=(7,5))
747 labels = 'miRNA RefSeq','Template', 'Unmapped','Non-template'
748 sizes = [c_mat_counts, c_tem_counts, c_unmap_counts,c_non_counts]
749 colors = ['gold', 'yellowgreen', 'lightcoral', 'lightskyblue']
750 ax1 = plt.subplot2grid((1,2),(0,0))
751 patches, texts, autotexts=plt.pie(sizes, labels=labels, colors=colors, startangle=140,autopct='%1.1f%%',radius=0.8)
752 [x.set_fontsize(8) for x in texts]
753 plt.title('Control Group (reads)',fontsize=12)
754 labels = 'miRNA RefSeq','Template', 'Unmapped','non-template'
755 sizes = [t_mat_counts, t_tem_counts, t_unmap_counts, t_non_counts]
756 colors = ['gold', 'yellowgreen', 'lightcoral', 'lightskyblue']
757 ax2 = plt.subplot2grid((1,2),(0,1))
758 patches, texts, autotexts=plt.pie(sizes, labels=labels, colors=colors, startangle=140,autopct='%1.1f%%',radius=0.8)
759 [x.set_fontsize(8) for x in texts]
760 plt.title('Treated Group (reads)', fontsize=12)
761 plt.savefig('pie_non.png',dpi=300)
762
763 ######################################################################################################################################################
764
765 def merging_names(LH2E_copy,new):
766
767 dupes=[]
768 final_LH2E =[]
769
770 for num in range(len(LH2E_copy)):
771
772 if LH2E_copy[num][1] not in final_LH2E and LH2E_copy[num][0] not in final_LH2E:
773 final_LH2E.append(LH2E_copy[num][1])
774 final_LH2E.append(LH2E_copy[num][0])
775 else:
776 dupes.append(LH2E_copy[num][1])
777
778 dupes=list(set(dupes))
779
780 for i in range(len(dupes)):
781 dupes[i]=[dupes[i]]
782
783 for x in LH2E_copy:
784 for y in dupes:
785 if x[1]==y[0]:
786 fl=0
787 if len(y)==1:
788 y.append(x[0])
789 else:
790 for i in range(1,len(y)):
791 if y[i].split("_")[0]==x[0].split("_")[0]:
792 fl=1
793 if len(x[0])<len(y[i]):
794 del y[i]
795 y.append(x[0])
796 break
797
798 if fl==0:
799 y.append((x[0]))
800
801 for y in dupes:
802 if len(y)>2:
803 for i in range(len(y)-1,1,-1):
804 y[1]=y[1]+"/"+y[i]
805 del y[i]
806
807
808 for x in LH2E_copy:
809 for y in dupes:
810 if x[1]==y[0]:
811 x[0]=y[1]
812
813 LH2E_copy.sort()
814 LH2E_copy=list(LH2E_copy for LH2E_copy,_ in itertools.groupby(LH2E_copy))
815
816 new.extend(LH2E_copy)
817
818
819 ######################################################################################################################################################
820 def pie_temp(merge_LH2E,c_unmap,c_unmap_counts,merge_LH8E,t_unmap,t_unmap_counts):
821
822 c_samples=[[x[0],x[1],sum(int(i) for i in x[2:])] for x in merge_LH2E]
823 t_samples=[[x[0],x[1],sum(int(i) for i in x[2:])] for x in merge_LH8E]
824
825 c_templ = 0
826 c_tem_counts = 0
827 c_mature = 0
828 c_mat_counts = 0
829 t_templ = 0
830 t_tem_counts = 0
831 t_mature = 0
832 t_mat_counts = 0
833
834 for x in c_samples:
835
836 if "/" not in x[0]:
837 if "chr" in x[0].split("_")[-1]:
838 c_mature+=1
839 c_mat_counts += x[2]
840 else:
841 c_templ+=1
842 c_tem_counts += x[2]
843 else:
844 f=0
845 for y in x[0].split("/"):
846 if "chr" in y.split("_")[-1]:
847 c_mature+=1
848 c_mat_counts += x[2]
849 f=1
850 break
851 if f==0:
852 c_templ+=1
853 c_tem_counts += x[2]
854
855 for x in t_samples:
856
857 if "/" not in x[0]:
858 if "chr" in x[0].split("_")[-1]:
859 t_mature+=1
860 t_mat_counts += x[2]
861 else:
862 t_templ+=1
863 t_tem_counts += x[2]
864 else:
865 f=0
866 for y in x[0].split("/"):
867 if "chr" in y.split("_")[-1]:
868 t_mature+=1
869 t_mat_counts += x[2]
870 f=1
871 break
872 if f==0:
873 t_templ+=1
874 t_tem_counts += x[2]
875
876
877 fig = plt.figure()
878 labels = 'miRNA RefSeq','Template', 'Unmapped'
879 sizes = [c_mat_counts, c_tem_counts, c_unmap_counts]
880 colors = ['gold', 'yellowgreen', 'lightskyblue']
881 explode = (0.2, 0.05, 0.1)
882 ax1 = plt.subplot2grid((1,2),(0,0))
883 patches, texts, autotexts=plt.pie(sizes, labels=labels, colors=colors, startangle=140,autopct='%1.1f%%',radius=0.8)
884 [x.set_fontsize(8) for x in texts]
885 plt.title('Control group (reads)', fontsize=12)
886 labels = 'miRNA RefSeq','Template', 'Unmapped'
887 sizes = [t_mat_counts, t_tem_counts, t_unmap_counts]
888 colors = ['gold', 'yellowgreen', 'lightskyblue']
889 explode = (0.2, 0.05, 0.1)
890 ax2 = plt.subplot2grid((1,2),(0,1))
891 patches, texts, autotexts=plt.pie(sizes, labels=labels, colors=colors, startangle=140,autopct='%1.1f%%',radius=0.8)
892 [x.set_fontsize(8) for x in texts]
893 plt.title('Treated group (reads)',fontsize = 12)
894 plt.savefig('pie_tem.png',dpi=300)
895
896 ###################################################################################################################################################################################################################
897
898 def make_spider(merge_LH2E,merge_LH8E):
899
900 c_samples=[[x[0],x[1],sum(int(i) for i in x[2:])] for x in merge_LH2E]
901 t_samples=[[x[0],x[1],sum(int(i) for i in x[2:])] for x in merge_LH8E]
902
903 c_5 = 0
904 c_5_counts = 0
905 c_3 = 0
906 c_3_counts = 0
907 c_both =0
908 c_both_counts=0
909 c_mature = 0
910 c_mat_counts = 0
911 c_exception=0
912 c_exception_counts=0
913
914
915 t_5 = 0
916 t_5_counts = 0
917 t_3 = 0
918 t_3_counts = 0
919 t_both = 0
920 t_both_counts = 0
921 t_mature = 0
922 t_mat_counts = 0
923 t_exception = 0
924 t_exception_counts=0
925
926 for x in c_samples:
927
928 if "/" not in x[0]:
929 if "chr" in x[0].split("_")[-1]:
930 c_mature+=1
931 c_mat_counts += x[2]
932 elif 0 == int(x[0].split("_")[-1]):
933 c_5+=1
934 c_5_counts += x[2]
935 elif 0 == int(x[0].split("_")[-2]):
936 c_3+=1
937 c_3_counts += x[2]
938 else:
939 c_both+=1
940 c_both_counts+=x[2]
941
942 else:
943 f=0
944 for y in x[0].split("/"):
945 if "chr" in y.split("_")[-1]:
946 c_mature+=1
947 c_mat_counts += x[2]
948 f=1
949 break
950 if f==0:
951 for y in x[0].split("/"):
952 c_exception+=1
953 c_exception_counts += x[2]
954
955
956 for x in t_samples:
957
958 if "/" not in x[0]:
959 if "chr" in x[0].split("_")[-1]:
960 t_mature+=1
961 t_mat_counts += x[2]
962 elif 0 == int(x[0].split("_")[-1]):
963 t_5+=1
964 t_5_counts += x[2]
965 elif 0 == int(x[0].split("_")[-2]):
966 t_3+=1
967 t_3_counts += x[2]
968 else:
969 t_both+=1
970 t_both_counts+=x[2]
971
972 else:
973 f=0
974 for y in x[0].split("/"):
975 if "chr" in y.split("_")[-1]:
976 t_mature+=1
977 t_mat_counts += x[2]
978 f=1
979 break
980 if f==0:
981 for y in x[0].split("/"):
982 t_exception+=1
983 t_exception_counts += x[2]
984
985
986 c_all = c_5+c_3+c_both+c_mature+c_exception
987 c_all_counts = c_5_counts + c_3_counts + c_both_counts + c_mat_counts + c_exception_counts
988
989 t_all = t_5+t_3+t_both+t_mature + t_exception
990 t_all_counts = t_5_counts + t_3_counts + t_both_counts + t_mat_counts + t_exception_counts
991
992 c_5 = round(c_5/c_all*100,2)
993 c_3 = round(c_3/c_all*100,2)
994 c_both = round(c_both/c_all*100,2)
995 c_mature = round(c_mature/c_all*100,2)
996 c_exception = round(c_exception/c_all*100,2)
997
998 c_5_counts = round(c_5_counts/c_all_counts*100,2)
999 c_3_counts = round(c_3_counts/c_all_counts*100,2)
1000 c_both_counts = round(c_both_counts/c_all_counts*100,2)
1001 c_mat_counts = round(c_mat_counts/c_all_counts*100,2)
1002 c_exception_counts = round(c_exception_counts/c_all_counts*100,2)
1003
1004 t_5 = round(t_5/t_all*100,2)
1005 t_3 = round(t_3/t_all*100,2)
1006 t_both = round(t_both/t_all*100,2)
1007 t_mature = round(t_mature/t_all*100,2)
1008 t_exception = round(t_exception/t_all*100,2)
1009
1010 t_5_counts = round(t_5_counts/t_all_counts*100,2)
1011 t_3_counts = round(t_3_counts/t_all_counts*100,2)
1012 t_both_counts = round(t_both_counts/t_all_counts*100,2)
1013 t_mat_counts = round(t_mat_counts/t_all_counts*100,2)
1014 t_exception_counts = round(t_exception_counts/t_all_counts*100,2)
1015
1016 radar_max = max(c_5, c_3, c_both,c_mature,c_exception,t_5,t_3,t_both,t_mature,t_exception)
1017 radar_max_counts = max(c_5_counts,c_3_counts,c_both_counts,c_mat_counts,c_exception_counts,t_5_counts,t_3_counts,t_both_counts,t_mat_counts,t_exception_counts)
1018
1019 df=pd.DataFrame({
1020 'group':['Controls','Treated'],
1021 """5' and 3' isomiRs""":[c_both,t_both],
1022 """3' isomiRs""":[c_3,t_3],
1023 'miRNA RefSeq':[c_mature,t_mature],
1024 """5' isomiRs""":[c_5,t_5],
1025 'Others*':[c_exception,t_exception]})
1026
1027 df1=pd.DataFrame({
1028 'group':['Controls','Treated'],
1029 """5' and 3' isomiRs""":[c_both_counts,t_both_counts],
1030 """3' isomiRs""":[c_3_counts,t_3_counts],
1031 'miRNA RefSeq':[c_mat_counts,t_mat_counts],
1032 """5' isomiRs""":[c_5_counts,t_5_counts],
1033 'Others*':[c_exception_counts,t_exception_counts]})
1034
1035 spider_last(df,radar_max,1)
1036 spider_last(df1,radar_max_counts,2)
1037
1038
1039 def spider_last(df,radar_max,flag):
1040 # ------- PART 1: Create background
1041 fig = plt.figure()
1042 # number of variable
1043 categories=list(df)[1:]
1044 N = len(categories)
1045
1046 # What will be the angle of each axis in the plot? (we divide the plot / number of variable)
1047 angles = [n / float(N) * 2 * pi for n in range(N)]
1048 angles += angles[:1]
1049
1050 # Initialise the spider plot
1051 ax = plt.subplot(111, polar=True)
1052
1053 # If you want the first axis to be on top:
1054 ax.set_theta_offset(pi/2)
1055 ax.set_theta_direction(-1)
1056
1057 # Draw one axe per variable + add labels labels yet
1058 plt.xticks(angles[:-1], categories, fontsize=11)
1059
1060 # Draw ylabels
1061 radar_max=round(radar_max+radar_max*0.1)
1062 mul=len(str(radar_max))-1
1063 maxi=int(math.ceil(radar_max / pow(10,mul))) * pow(10,mul)
1064 sep = round(maxi/4)
1065 plt.yticks([sep, 2*sep, 3*sep, 4*sep, 5*sep], [str(sep)+'%', str(2*sep)+'%', str(3*sep)+'%', str(4*sep)+'%', str(5*sep)+'%'], color="grey", size=10)
1066 plt.ylim(0, maxi)
1067
1068 # ------- PART 2: Add plots
1069
1070 # Plot each individual = each line of the data
1071 # I don't do a loop, because plotting more than 3 groups makes the chart unreadable
1072
1073 # Ind1
1074 values=df.loc[0].drop('group').values.flatten().tolist()
1075 values += values[:1]
1076 ax.plot(angles, values,'-o', linewidth=1, linestyle='solid', label="Controls")
1077 ax.fill(angles, values, 'b', alpha=0.1)
1078
1079 # Ind2
1080 values=df.loc[1].drop('group').values.flatten().tolist()
1081 values += values[:1]
1082 ax.plot(angles, values, '-o' ,linewidth=1, linestyle='solid', label="Treated")
1083 ax.fill(angles, values, 'r', alpha=0.1)
1084
1085 # Add legend
1086 if flag==1:
1087 plt.legend(loc='upper right', bbox_to_anchor=(0.0, 0.1))
1088 plt.savefig('spider_non_red.png',dpi=300)
1089 else:
1090 plt.legend(loc='upper right', bbox_to_anchor=(0.0, 0.1))
1091 plt.savefig('spider_red.png',dpi=300)
1092
1093
1094 #############################################################################################################################################################################################################
1095 def hist_red(samples,flag):
1096 lengths=[]
1097 cat=[]
1098 total_reads=0
1099 seq=[]
1100
1101 if flag == "c":
1102 title = "Length Distribution of Control group (Redudant reads)"
1103 if flag == "t":
1104 title = "Length Distribution of Treated group (Redudant reads)"
1105
1106 for i in samples:
1107 for x in i:
1108 lengths.append(len(x[9]))
1109 if x[1]=="0":
1110 seq.append([x[9],x[0].split("-")[1],"Mapped"])
1111 cat.append("Mapped")
1112 if x[1] == "4":
1113 seq.append([x[9],x[0].split("-")[1],"Unmapped"])
1114 cat.append("Unmapped")
1115
1116 uni_len=list(set(lengths))
1117 uni_len=[x for x in uni_len if x<=35]
1118 low=min(lengths)
1119 up=max(lengths)
1120 seq.sort()
1121 uni_seq=list(seq for seq,_ in itertools.groupby(seq))
1122 dim=up-low
1123
1124 if dim>20:
1125 s=5
1126 else:
1127 s=8
1128
1129 total_reads+=sum([int(x[1]) for x in uni_seq])
1130
1131 map_reads=[]
1132 unmap_reads=[]
1133 length=[]
1134 for y in uni_len:
1135 map_temp=0
1136 unmap_temp=0
1137 for x in uni_seq:
1138 if len(x[0])==y and x[2]=="Mapped":
1139 map_temp+=int(x[1])
1140 if len(x[0])==y and x[2]=="Unmapped":
1141 unmap_temp+=int(x[1])
1142 if y<=35:
1143 length.append(y)
1144 map_reads.append(round(map_temp/total_reads*100,2))
1145 unmap_reads.append(round(unmap_temp/total_reads*100,2))
1146
1147 ylim=max([sum(x) for x in zip(unmap_reads, map_reads)])
1148 ylim=ylim+ylim*20/100
1149 fig, ax = plt.subplots()
1150 width=0.8
1151 ax.bar(length, unmap_reads, width, label='Unmapped')
1152 h=ax.bar(length, map_reads, width, bottom=unmap_reads, label='Mapped')
1153 plt.xticks(np.arange(length[0], length[-1]+1, 1))
1154 plt.xlabel('Length (nt)',fontsize=14)
1155 plt.ylabel('Percentage',fontsize=14)
1156 plt.title(title,fontsize=14)
1157 ax.legend()
1158 plt.ylim([0, ylim])
1159 ax.grid(axis='y',linewidth=0.2)
1160
1161 if flag=='c':
1162 plt.savefig('c_hist_red.png',dpi=300)
1163
1164 if flag=='t':
1165 plt.savefig('t_hist_red.png',dpi=300)
1166
1167 #################################################################################################################
1168
1169
1170 def logo_seq_red(merge, flag):
1171
1172 if flag=="c":
1173 titlos="Control group (Redundant)"
1174 file_logo="c_logo.png"
1175 file_bar="c_bar.png"
1176 if flag=="t":
1177 titlos="Treated group (Redundant)"
1178 file_logo="t_logo.png"
1179 file_bar="t_bar.png"
1180
1181 c_samples=[[x[0],x[1],sum(int(i) for i in x[2:])] for x in merge]
1182
1183 A=[0]*5
1184 C=[0]*5
1185 G=[0]*5
1186 T=[0]*5
1187 total_reads=0
1188
1189 for y in c_samples:
1190 if "/" in y[0]:
1191 length=[]
1192 for x in y[0].split("/"):
1193 length.append([len(x.split("_")[-1]),x.split("_")[-1],y[2]])
1194
1195 best=min(length)
1196 total_reads+=best[2]
1197 for i in range(5):
1198 if i<len(best[1]):
1199 if best[1][i] == "A":
1200 A[i]+=best[2]
1201 elif best[1][i] == "C":
1202 C[i]+=best[2]
1203 elif best[1][i] == "G":
1204 G[i]+=best[2]
1205 else:
1206 T[i]+=best[2]
1207 else:
1208 total_reads+=y[2]
1209 for i in range(5):
1210 if i<len(y[0].split("_")[-1]):
1211 if y[0].split("_")[-1][i] == "A":
1212 A[i]+=(y[2])
1213 elif y[0].split("_")[-1][i] == "C":
1214 C[i]+=(y[2])
1215 elif y[0].split("_")[-1][i] == "G":
1216 G[i]+=(y[2])
1217 else:
1218 T[i]+=y[2]
1219
1220 A[:] = [round(x*100,1) / total_reads for x in A]
1221 C[:] = [round(x*100,1) / total_reads for x in C]
1222 G[:] = [round(x*100,1) / total_reads for x in G]
1223 T[:] = [round(x*100,1) / total_reads for x in T]
1224
1225
1226
1227 data = {'A':A,'C':C,'G':G,'T':T}
1228 df = pd.DataFrame(data, index=[1,2,3,4,5])
1229 h=df.plot.bar(color=tuple(["g", "b","gold","r"]) )
1230 h.grid(axis='y',linewidth=0.2)
1231 plt.xticks(rotation=0, ha="right")
1232 plt.ylabel("Counts (%)",fontsize=18)
1233 plt.xlabel("Positions (nt)",fontsize=18)
1234 plt.title(titlos,fontsize=20)
1235 plt.tight_layout()
1236 plt.savefig(file_bar, dpi=300)
1237
1238 import logomaker as lm
1239 crp_logo = lm.Logo(df, font_name = 'DejaVu Sans')
1240 crp_logo.style_spines(visible=False)
1241 crp_logo.style_spines(spines=['left', 'bottom'], visible=True)
1242 crp_logo.style_xticks(rotation=0, fmt='%d', anchor=0)
1243
1244 # style using Axes methods
1245 crp_logo.ax.set_title(titlos,fontsize=18)
1246 crp_logo.ax.set_ylabel("Counts (%)", fontsize=16,labelpad=5)
1247 crp_logo.ax.set_xlabel("Positions (nt)",fontsize=16, labelpad=5)
1248 crp_logo.ax.xaxis.set_ticks_position('none')
1249 crp_logo.ax.xaxis.set_tick_params(pad=-1)
1250 figure = plt.gcf()
1251 figure.set_size_inches(6, 4)
1252 crp_logo.fig.savefig(file_logo,dpi=300)
1253
1254 ##########################################################################################################################################################################################################
1255
1256
1257
1258 def logo_seq_non_red(merge_LH2E):
1259
1260 c_samples=[[x[0],x[1],sum(int(i) for i in x[2:])] for x in merge_LH2E]
1261
1262 A=[0]*5
1263 C=[0]*5
1264 G=[0]*5
1265 T=[0]*5
1266
1267 for y in c_samples:
1268 if "/" in y[0]:
1269 length=[]
1270 for x in y[0].split("/"):
1271 length.append([len(x.split("_")[-1]),x.split("_")[-1],y[2]])
1272
1273 best=min(length)
1274 for i in range(5):
1275 if i<len(best[1]):
1276 if best[1][i] == "A":
1277 A[i]+=1
1278 elif best[1][i] == "C":
1279 C[i]+=1
1280 elif best[1][i] == "G":
1281 G[i]+=1
1282 else:
1283 T[i]+=1
1284 else:
1285 for i in range(5):
1286 if i<len(y[0].split("_")[-1]):
1287 if y[0].split("_")[-1][i] == "A":
1288 A[i]+=1
1289 elif y[0].split("_")[-1][i] == "C":
1290 C[i]+=1
1291 elif y[0].split("_")[-1][i] == "G":
1292 G[i]+=1
1293 else:
1294 T[i]+=1
1295
1296 data = {'A':A,'C':C,'G':G,'T':T}
1297 df = pd.DataFrame(data, index=[1,2,3,4,5])
1298 h=df.plot.bar(title="Non-templated nucleotides after templated sequence",color=tuple(["g", "b","gold","r"]))
1299 h.set_xlabel("Positions (nt)")
1300 h.set_ylabel("Unique sequences")
1301 plt.xticks(rotation=0, ha="right")
1302 plt.tight_layout()
1303 plt.savefig("bar2.png", dpi=300)
1304
1305
1306 import logomaker as lm
1307 crp_logo = lm.Logo(df, font_name = 'DejaVu Sans')
1308
1309 # style using Logo methods
1310 crp_logo.style_spines(visible=False)
1311 crp_logo.style_spines(spines=['left', 'bottom'], visible=True)
1312 crp_logo.style_xticks(rotation=0, fmt='%d', anchor=0)
1313
1314 # style using Axes methods
1315 crp_logo.ax.set_ylabel("Unique sequences", labelpad=5)
1316 crp_logo.ax.set_xlabel("Positions (nt)", labelpad=5)
1317 crp_logo.ax.xaxis.set_ticks_position('none')
1318 crp_logo.ax.xaxis.set_tick_params(pad=-1)
1319 crp_logo.ax.set_title("Non-redundant")
1320 figure = plt.gcf()
1321 crp_logo.fig.savefig('logo2.png', dpi=300)
1322
1323
1324 ###################################################################################################################################################################################################################
1325
1326 def ssamples1(tem_names,tem_samp,non_names,non_samp,folder,pro):
1327
1328 for i in range(2,len(tem_samp[0])):
1329
1330 fp = open(folder+tem_names[i-2]+'.txt','w')
1331 fp.write("miRNA id"+"\t"+tem_names[i-2]+"\n")
1332
1333 for x in tem_samp:
1334 fp.write("%s" % "\t".join([x[0],x[i]])+"\n")
1335
1336 for j in range(len(non_names)):
1337 if non_names[j]==tem_names[i-2]:
1338 for x in non_samp:
1339 fp.write("%s" % "\t".join([x[0],x[j+2]])+"\n")
1340 fp.close()
1341
1342 ###################################################################################################################################################################################################################
1343
1344 def download_matures(matures,org_name):
1345
1346 #url = 'ftp://mirbase.org/pub/mirbase/21/mature.fa.gz'
1347 url = 'ftp://mirbase.org/pub/mirbase/CURRENT/mature.fa.gz'
1348 data = urllib.request.urlopen(url).read()
1349 file_mirna = gzip.decompress(data).decode('utf-8')
1350 file_mirna = file_mirna.split("\n")
1351
1352 for i in range(0,len(file_mirna)-1,2):
1353
1354 if org_name in file_mirna[i]:
1355 matures.append(file_mirna[i])
1356 matures.append(file_mirna[i+1])
1357
1358 ###################################################################################################################################################################################################################
1359 def non_template_ref(sc,st,all_isoforms):
1360
1361 pre_uni_seq_con = list(sc)
1362 pre_uni_seq_tre = list(st)
1363
1364 for x in pre_uni_seq_con:
1365 for y in x:
1366 if ">"+y[2] not in all_isoforms and ")_" in y[2] :
1367 all_isoforms.append(">"+y[2])
1368 all_isoforms.append(y[9])
1369
1370
1371 for x in pre_uni_seq_tre:
1372 for y in x:
1373 if ">"+y[2] not in all_isoforms and ")_" in y[2]:
1374 all_isoforms.append(">"+y[2])
1375 all_isoforms.append(y[9])
1376
1377 ################################################################################################################################################################################################
1378
1379 def deseqe2(sample,mir_names,l,new_d,sample_name,sample_order):
1380
1381 for y in mir_names:
1382 flag=0
1383 for x in sample:
1384 if y[0]==x[0]:
1385 flag=1
1386 break
1387 if flag==0:
1388 sample.append([y[0],"0",y[1]])
1389
1390 sample.sort(key=lambda x: x[0])
1391 sample=list(sample for sample,_ in itertools.groupby(sample))
1392
1393 l.acquire()
1394 new_d.append(sample)
1395 sample_order.append(sample_name)
1396 l.release()
1397
1398 ###############################################################################################################################################################################################
1399
1400 if __name__ == '__main__':
1401
1402 starttime = time.time()
1403
1404 q1 = Queue()
1405 q2 = Queue()
1406 lock = Lock()
1407 manager = Manager()
1408
1409 mature_mirnas=manager.list()
1410 ps_mature=Process(target=download_matures,args=(mature_mirnas,args.org_name))
1411 ps_mature.start()
1412
1413 args.control[0]=args.control[0][1:]
1414 args.control[len(args.control)-1][:-1]
1415 control = [(args.control[i:i+2]) for i in range(0, len(args.control), 2)]
1416
1417 args.treated[0]=args.treated[0][1:]
1418 args.treated[len(args.treated)-1][:-1]
1419 treated = [(args.treated[i:i+2]) for i in range(0, len(args.treated), 2)]
1420
1421
1422 ############## Detection of templated isoforms ################
1423
1424 radar = manager.list([0,0,0,0])
1425 samples = manager.list()
1426 data= manager.list()
1427 names_con=manager.list()
1428 samples_mirna_names=manager.list()
1429 deseq=manager.list()
1430 unmap_seq=manager.Value('i',0)
1431 unmap_counts=manager.Value('i',0)
1432 LH2E_names=manager.list()
1433 ini_c_samples = manager.list()
1434
1435
1436 radar1 = manager.list([0,0,0,0])
1437 samples1 = manager.list()
1438 data1 = manager.list()
1439 names_tre = manager.list()
1440 samples_mirna_names1=manager.list()
1441 deseq1=manager.list()
1442 unmap1_seq = manager.Value('i',0)
1443 unmap1_counts = manager.Value('i',0)
1444 LH8E_names=manager.list()
1445 ini_t_samples = manager.list()
1446 ps_mature.join()
1447
1448
1449 mature_mirnas=list(mature_mirnas)
1450
1451
1452 starttime1 = time.time()
1453 ps_sam = [Process(target=sam,args=(mature_mirnas,path[1][:-1],path[0].split(",")[0],"c",lock,samples,data,names_con,unmap_seq,samples_mirna_names,deseq,LH2E_names,"0",ini_c_samples,unmap_counts)) for path in control]
1454 ps_sam.extend([Process(target=sam,args=(mature_mirnas,path[1][:-1],path[0].split(",")[0],"t",lock,samples1,data1,names_tre,unmap1_seq,samples_mirna_names1,deseq1,LH8E_names,"0",ini_t_samples,unmap1_counts)) for path in treated])
1455
1456 [p.start() for p in ps_sam]
1457 [p.join() for p in ps_sam]
1458 print('SAM took {} seconds'.format(time.time() - starttime1))
1459
1460 ps_hist=[Process(target=hist_red,args=(ini_c_samples,'c'))]
1461 ps_hist.extend([Process(target=hist_red,args=(ini_t_samples,'t'))])
1462 [x.start() for x in ps_hist]
1463
1464 starttime200=time.time()
1465
1466 sc = list(samples)
1467 st = list(samples1)
1468
1469 names_con=list(names_con)
1470 names_tre=list(names_tre)
1471 samples_mirna_names=list(samples_mirna_names)
1472 samples_mirna_names.sort()
1473 samples_mirna_names=list(samples_mirna_names for samples_mirna_names,_ in itertools.groupby(samples_mirna_names))
1474
1475 samples_mirna_names1=list(samples_mirna_names1)
1476 samples_mirna_names1.sort()
1477 samples_mirna_names1=list(samples_mirna_names1 for samples_mirna_names1,_ in itertools.groupby(samples_mirna_names1))
1478
1479 deseq=list(deseq)
1480 deseq1=list(deseq1)
1481
1482 new_names_con=manager.list()
1483 new_names_tre=manager.list()
1484 new_deseq=manager.list()
1485 new_deseq1=manager.list()
1486 ps_deseq=[Process(target=deseqe2,args=(sampp,samples_mirna_names,lock,new_deseq,names_con[i],new_names_con)) for i,sampp in enumerate(deseq)]
1487 ps_deseq.extend([Process(target=deseqe2,args=(sampp,samples_mirna_names1,lock,new_deseq1,names_tre[i],new_names_tre)) for i,sampp in enumerate(deseq1)])
1488
1489 [z.start() for z in ps_deseq]
1490 [z.join() for z in ps_deseq]
1491 new_deseq=list(new_deseq)
1492 new_deseq1=list(new_deseq1)
1493
1494 LH2E=[[x[0],x[2]] for x in new_deseq[0]]
1495 [LH2E[i].append(y[i][1]) for i,_ in enumerate(LH2E) for y in new_deseq]
1496
1497 LH8E=[[x[0],x[2]] for x in new_deseq1[0]]
1498 [LH8E[i].append(y[i][1]) for i,_ in enumerate(LH8E) for y in new_deseq1]
1499
1500 print('Deseq took {} seconds'.format(time.time() - starttime200))
1501
1502 merg_nam_LH2E=manager.list()
1503 merg_nam_LH8E=manager.list()
1504
1505 LH2E_copy=copy.deepcopy(list(LH2E))
1506 LH8E_copy=copy.deepcopy(list(LH8E))
1507
1508 fil_sort_tre=manager.list()
1509 fil_sort_con=manager.list()
1510 raw_sort_tre=manager.list()
1511 raw_sort_con=manager.list()
1512
1513 ps_main = Process(target=main_temp,args=(list(LH2E), samples_mirna_names, list(LH8E), samples_mirna_names1,1,list(names_con),list(names_tre),fil_sort_tre,fil_sort_con,raw_sort_tre,raw_sort_con))
1514 ps_main.start()
1515
1516 if args.anal=="2":
1517 all_iso = manager.list()
1518 ps_non_iso = Process(target=non_template_ref,args=(sc,st,all_iso))
1519 ps_non_iso.start()
1520
1521 ps_merge = [Process(target=merging_names,args=(LH2E_copy,merg_nam_LH2E))]
1522 ps_merge.extend([Process(target=merging_names,args=(LH8E_copy,merg_nam_LH8E))])
1523 [x.start() for x in ps_merge]
1524 [x.join() for x in ps_merge]
1525
1526 merg_nam_LH2E=list(merg_nam_LH2E)
1527 merg_nam_LH8E=list(merg_nam_LH8E)
1528
1529 starttime2 = time.time()
1530 procs = [Process(target=DB_write,args=(x[0],x[1],x[2],x[3],1)) for x in data]
1531 procs.extend([Process(target=DB_write,args=(x[0],x[1],x[2],x[3],1)) for x in data1])
1532 procs.extend([Process(target=make_spider,args=(merg_nam_LH2E,merg_nam_LH8E))])
1533 if args.anal == "1":
1534 procs.extend([Process(target=pie_temp,args=(merg_nam_LH2E,unmap_seq.value,unmap_counts.value,merg_nam_LH8E,unmap1_seq.value,unmap1_counts.value))])
1535
1536 [p.start() for p in procs]
1537
1538
1539 if args.anal=="1":
1540 [x.join() for x in ps_hist]
1541 [p.join() for p in procs]
1542 ps_pdf = Process(target=pdf_before_DE,args=(args.anal))
1543 ps_pdf.start()
1544
1545 print('Graphs took {} seconds'.format(time.time() - starttime2))
1546
1547 ps_main.join()
1548
1549 fil_sort_con=list(fil_sort_con)
1550 fil_sort_tre=list(fil_sort_tre)
1551 if fil_sort_con==[]:
1552 fil_sort_con=raw_sort_con
1553 fil_sort_tre=raw_sort_tre
1554
1555 raw_sort_con=list(raw_sort_con)
1556 raw_sort_tre=list(raw_sort_tre)
1557 names_con=list(new_names_con)
1558 names_tre=list(new_names_tre)
1559
1560 ps_write = Process(target=write_main,args=(raw_sort_con, raw_sort_tre, fil_sort_con, fil_sort_tre, names_con,names_tre,1))
1561 ps_write.start()
1562
1563 ps1_matrix = [Process(target=ssamples,args=(names_con,fil_sort_con,"Diff/temp_con/",0))]
1564 ps1_matrix.extend([Process(target=ssamples,args=(names_tre,fil_sort_tre,"Diff/temp_tre/",0))])
1565 [p.start() for p in ps1_matrix]
1566
1567 if args.anal=="1":
1568 ps_pdf.join()
1569 if args.anal=="2":
1570 [p.join() for p in procs]
1571 [x.join() for x in ps_hist]
1572
1573 ps_write.join()
1574 [p.join() for p in ps1_matrix]
1575
1576
1577 ############################## Detection of Both #######################################
1578
1579 starttime10 = time.time()
1580
1581 if args.anal == "2":
1582
1583 n_data= manager.list()
1584 n_names_con=manager.list()
1585 n_samples_mirna_names=manager.list()
1586 n_deseq=manager.list()
1587 n_LH2E_names=manager.list()
1588
1589 n_data1 = manager.list()
1590 n_names_tre = manager.list()
1591 n_samples_mirna_names1=manager.list()
1592 n_deseq1=manager.list()
1593 n_LH8E_names=manager.list()
1594
1595
1596 new_mat_mirnas = list(mature_mirnas)
1597 ps_non_iso.join()
1598
1599 all_iso=list(all_iso)
1600 new_mat_mirnas.extend(all_iso)
1601
1602 starttime11=time.time()
1603
1604 ps_sam = [Process(target=non_sam,args=(new_mat_mirnas,path[1][:-1],path[0].split(",")[0],"c",lock,n_data,n_names_con,n_deseq,n_samples_mirna_names,n_LH2E_names)) for path in control]
1605 ps_sam.extend([Process(target=non_sam,args=(new_mat_mirnas,path[1][:-1],path[0].split(",")[0],"t",lock,n_data1,n_names_tre,n_deseq1,n_samples_mirna_names1,n_LH8E_names)) for path in treated])
1606
1607 [p.start() for p in ps_sam]
1608 [p.join() for p in ps_sam]
1609
1610 print('Non-sam took {} seconds'.format(time.time() - starttime11))
1611
1612 starttime12=time.time()
1613
1614 n_names_con=list(n_names_con)
1615 n_names_tre=list(n_names_tre)
1616 n_samples_mirna_names=list(n_samples_mirna_names)
1617 n_samples_mirna_names.sort()
1618 n_samples_mirna_names=list(n_samples_mirna_names for n_samples_mirna_names,_ in itertools.groupby(n_samples_mirna_names))
1619
1620 n_samples_mirna_names1=list(n_samples_mirna_names1)
1621 n_samples_mirna_names1.sort()
1622 n_samples_mirna_names1=list(n_samples_mirna_names1 for n_samples_mirna_names1,_ in itertools.groupby(n_samples_mirna_names1))
1623
1624 n_deseq=list(n_deseq)
1625 n_deseq1=list(n_deseq1)
1626
1627 new_n_names_con=manager.list()
1628 new_n_names_tre=manager.list()
1629 n_new_deseq=manager.list()
1630 n_new_deseq1=manager.list()
1631 ps_deseq=[Process(target=deseqe2,args=(sampp,n_samples_mirna_names,lock,n_new_deseq,n_names_con[i],new_n_names_con)) for i,sampp in enumerate(n_deseq)]
1632 ps_deseq.extend([Process(target=deseqe2,args=(sampp,n_samples_mirna_names1,lock,n_new_deseq1,n_names_tre[i],new_n_names_tre)) for i,sampp in enumerate(n_deseq1)])
1633
1634 [x.start() for x in ps_deseq]
1635 [x.join() for x in ps_deseq]
1636 n_new_deseq=list(n_new_deseq)
1637 n_new_deseq1=list(n_new_deseq1)
1638
1639 print([len(n_new_deseq[0]),len(n_new_deseq[1])])
1640 print([len(n_new_deseq1[0]),len(n_new_deseq1[1])])
1641
1642 n_LH2E=[[x[0],x[2]] for x in n_new_deseq[0]]
1643 [n_LH2E[i].append(y[i][1]) for i,_ in enumerate(n_LH2E) for y in n_new_deseq]
1644
1645 n_LH8E=[[x[0],x[2]] for x in n_new_deseq1[0]]
1646 [n_LH8E[i].append(y[i][1]) for i,_ in enumerate(n_LH8E) for y in n_new_deseq1]
1647
1648 print('Non-deseq took {} seconds'.format(time.time() - starttime12))
1649
1650 merg_nam_n_LH2E=manager.list()
1651 merg_nam_n_LH8E=manager.list()
1652
1653 n_LH2E_copy=copy.deepcopy(list(n_LH2E))
1654 n_LH8E_copy=copy.deepcopy(list(n_LH8E))
1655
1656 n_sort_tre=manager.list()
1657 n_sort_con=manager.list()
1658
1659 n_fil_sort_con=manager.list()
1660 n_fil_sort_tre=manager.list()
1661 n_raw_sort_con=manager.list()
1662 n_raw_sort_tre=manager.list()
1663
1664 ps_main = Process(target=main_temp,args=(list(n_LH2E), n_samples_mirna_names, list(n_LH8E), n_samples_mirna_names1,1,list(n_names_con),list(n_names_tre),n_fil_sort_tre,n_fil_sort_con,n_raw_sort_tre,n_raw_sort_con))
1665 ps_main.start()
1666
1667 ps_merge = [Process(target=merging_names,args=(n_LH2E_copy,merg_nam_n_LH2E))]
1668 ps_merge.extend([Process(target=merging_names,args=(n_LH8E_copy,merg_nam_n_LH8E))])
1669 [p.start() for p in ps_merge]
1670 [p.join() for p in ps_merge]
1671
1672 merg_nam_n_LH2E=list(merg_nam_n_LH2E)
1673 merg_nam_n_LH8E=list(merg_nam_n_LH8E)
1674
1675 procs = [Process(target=DB_write,args=(x[0],x[1],x[2],x[3],2)) for x in n_data]
1676 procs.extend([Process(target=DB_write,args=(x[0],x[1],x[2],x[3],2)) for x in n_data1])
1677 procs.extend([Process(target=logo_seq_red,args=(merg_nam_n_LH2E,'c'))])
1678 procs.extend([Process(target=logo_seq_red,args=(merg_nam_n_LH8E,'t'))])
1679 procs.extend([Process(target=pie_non_temp,args=(merg_nam_LH2E,merg_nam_n_LH2E,merg_nam_LH8E,merg_nam_n_LH8E,unmap_seq.value,unmap1_seq.value,unmap_counts.value,unmap1_counts.value))])
1680
1681 starttime13=time.time()
1682 [p.start() for p in procs]
1683 [p.join() for p in procs]
1684
1685 print('Graphs took {} seconds'.format(time.time() - starttime13))
1686
1687 procs1 = Process(target=pdf_before_DE,args=(args.anal))
1688 procs1.start()
1689
1690 starttime14=time.time()
1691 ps_main.join()
1692
1693 n_fil_sort_con=list(n_fil_sort_con)
1694 n_fil_sort_tre=list(n_fil_sort_tre)
1695 if n_fil_sort_con==[]:
1696 n_fil_sort_con=n_raw_sort_con
1697 n_fil_sort_tre=n_raw_sort_tre
1698
1699 n_raw_sort_con=list(n_raw_sort_con)
1700 n_raw_sort_tre=list(n_raw_sort_tre)
1701 n_names_con=list(new_n_names_con)
1702 n_names_tre=list(new_n_names_tre)
1703
1704 ps_write = Process(target=write_main,args=(n_raw_sort_con, n_raw_sort_tre,n_fil_sort_con, n_fil_sort_tre, n_names_con, n_names_tre,2))
1705 ps_write.start()
1706
1707 ps1_matrix = [Process(target=ssamples1,args=(n_names_con,n_fil_sort_con,names_con,fil_sort_con,"Diff/n_temp_con/",0))]
1708 ps1_matrix.extend([Process(target=ssamples1,args=(n_names_tre,n_fil_sort_tre,names_tre,fil_sort_tre,"Diff/n_temp_tre/",0))])
1709 [p.start() for p in ps1_matrix]
1710
1711 ps_write.join()
1712 [p.join() for p in ps1_matrix]
1713 procs1.join()
1714 print('That took {} seconds'.format(time.time() - starttime10))
1715 print('That took {} seconds'.format(time.time() - starttime))
1716
1717
1718
1719
1720
1721
1722
1723