Mercurial > repos > glogobyte > isoread

comparison mirbase_ultra_v2.py @ 14:3ad9701c7749 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author glogobyte
date Thu, 22 Oct 2020 07:30:47 +0000
parents
children
comparison
equal deleted inserted replaced
13:f1b31ee97e91 14:3ad9701c7749
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 zeros=["0"]*(len(LH8E[0])-2)
399 [LH8E_add_names[i].extend(zeros) for i,_ in enumerate(LH8E_add_names)]
400 LH8E=LH8E+LH8E_add_names
401
402 zeros=["0"]*(len(LH2E[0])-2)
403 [LH2E_add_names[i].extend(zeros) for i,_ in enumerate(LH2E_add_names)]
404 LH2E=LH2E+LH2E_add_names
405
406 dupes=[]
407 final_LH2E =[]
408
409 for num,_ in enumerate(LH2E):
410
411 if LH2E[num][1] not in final_LH2E and LH2E[num][0] not in final_LH2E:
412 final_LH2E.append(LH2E[num][1])
413 final_LH2E.append(LH2E[num][0])
414 else:
415 dupes.append(LH2E[num][1])
416
417
418 dupes=list(set(dupes))
419
420 dupes=[[x] for x in dupes]
421
422 for x in LH2E:
423 for y in dupes:
424 if x[1]==y[0]:
425 fl=0
426 if len(y)==1:
427 y.append(x[0])
428 else:
429 for i in range(1,len(y)):
430 if y[i].split("_")[0]==x[0].split("_")[0]:
431 fl=1
432 if len(x[0])<len(y[i]):
433 del y[i]
434 y.append(x[0])
435 break
436
437 if fl==0:
438 y.append((x[0]))
439
440 for y in dupes:
441 if len(y)>2:
442 for i in range(len(y)-1,1,-1):
443 y[1]=y[1]+"/"+y[i]
444 del y[i]
445
446 for x in LH2E:
447 for y in dupes:
448 if x[1]==y[0]:
449 x[0]=y[1]
450
451 for x in LH8E:
452 for y in dupes:
453 if x[1]==y[0]:
454 x[0]=y[1]
455
456
457 LH2E.sort()
458 LH2E=list(LH2E for LH2E,_ in itertools.groupby(LH2E))
459
460 LH8E.sort()
461 LH8E=list(LH8E for LH8E,_ in itertools.groupby(LH8E))
462
463 LH8E_new=[]
464 LH2E_new=[]
465
466 if int(args.per)!=-1:
467 percent=int(args.per)/100
468
469 c_col_filter=round(percent*(len(LH2E[1])-2))
470 t_col_filter=round(percent*(len(LH8E[1])-2))
471
472 for i, _ in enumerate(LH2E):
473 c_cols=0
474 t_cols=0
475
476 c_cols=sum([1 for j in range(len(LH2E[i])-2) if int(LH2E[i][j+2])>=int(args.count)])
477 t_cols=sum([1 for j in range(len(LH8E[i])-2) if int(LH8E[i][j+2])>=int(args.count)])
478
479 if c_cols>=c_col_filter or t_cols>=t_col_filter:
480 LH8E_new.append(LH8E[i])
481 LH2E_new.append(LH2E[i])
482
483 filter_LH2E.extend(LH2E_new)
484 filter_LH8E.extend(LH8E_new)
485 raw_LH2E.extend(LH2E)
486 raw_LH8E.extend(LH8E)
487
488 ##################################################################################################################################################################################################################
489
490 def write_main(raw_LH2E, raw_LH8E, fil_LH2E, fil_LH8E, names_con, names_tre, flag):
491
492 if flag == 1 and int(args.per)!=-1:
493 fp = open('Counts/Filtered '+args.n2 +' Templated Counts', 'w')
494 fp.write("Name\t")
495 fp.write("Sequence")
496 for y in names_tre:
497 fp.write("\t"+y)
498
499 for x in fil_LH8E:
500 fp.write("\n%s" % "\t".join(x))
501 fp.close()
502
503 fp = open('Counts/Filtered '+args.n1+' Templated Counts', 'w')
504 fp.write("Name\t")
505 fp.write("Sequence")
506 for y in names_con:
507 fp.write("\t"+y)
508
509 for x in fil_LH2E:
510 fp.write("\n%s" % "\t".join(x))
511 fp.close()
512
513
514 if flag == 2 and int(args.per)!=-1:
515 fp = open('Counts/Filtered '+args.n2+' Non-Templated Counts', 'w')
516 fp.write("Name\t")
517 fp.write("Sequence")
518 for y in names_tre:
519 fp.write("\t"+y)
520
521
522 for x in fil_LH8E:
523 fp.write("\n%s" % "\t".join(x))
524 fp.close()
525
526 fp = open('Counts/Filtered '+args.n1+' Non-Templated Counts', 'w')
527 fp.write("Name\t")
528 fp.write("Sequence")
529 for y in names_con:
530 fp.write("\t"+y)
531
532 for x in fil_LH2E:
533 fp.write("\n%s" % "\t".join(x))
534 fp.close()
535
536
537 if flag == 1:
538 fp = open('Counts/Raw '+args.n2+' Templated Counts', 'w')
539 fp.write("Name\t")
540 fp.write("Sequence")
541 for y in names_tre:
542 fp.write("\t"+y)
543
544 for x in raw_LH8E:
545 fp.write("\n%s" % "\t".join(x))
546 fp.close()
547
548 fp = open('Counts/Raw '+args.n1+' Templated Counts', 'w')
549 fp.write("Name\t")
550 fp.write("Sequence")
551 for y in names_con:
552 fp.write("\t"+y)
553
554 for x in raw_LH2E:
555 fp.write("\n%s" % "\t".join(x))
556 fp.close()
557
558
559 if flag == 2:
560 fp = open('Counts/Raw '+args.n2+' Non-Templated Counts', 'w')
561 fp.write("Name\t")
562 fp.write("Sequence")
563 for y in names_tre:
564 fp.write("\t"+y)
565
566
567 for x in raw_LH8E:
568 fp.write("\n%s" % "\t".join(x))
569 fp.close()
570
571 fp = open('Counts/Raw '+args.n1+' Non-Templated Counts', 'w')
572 fp.write("Name\t")
573 fp.write("Sequence")
574 for y in names_con:
575 fp.write("\t"+y)
576
577 for x in raw_LH2E:
578 fp.write("\n%s" % "\t".join(x))
579 fp.close()
580
581
582 #########################################################################################################################################
583
584 def ssamples(names,samp,folder,pro):
585
586 for i in range(2,len(samp[0])):
587
588 fp = open(folder+names[i-2]+'.txt','w')
589 fp.write("miRNA id"+"\t"+names[i-2]+"\n")
590
591 for x in samp:
592 fp.write("%s" % "\t".join([x[0],x[i]])+"\n")
593 fp.close()
594
595 ##################################################################################################################
596
597 def DB_write(con,name,unique_seq,sorted_uni_arms,f):
598
599 if f==1:
600 # Write a txt file with all the information
601 if con=="c":
602 fp = open('split1/'+name, 'w')
603
604 fp.write("%s\t%-42s\t%s\n\n" % ("Number of Reads","Name of isomir","Sequence"))
605 if con=="t":
606 fp = open('split2/'+name, 'w')
607 fp.write("%s\t%-42s\t%s\n\n" % ("Number of Reads","Name of isomir","Sequence"))
608
609
610 for i in range(len(sorted_uni_arms)):
611 temp = []
612 for j in range(len(unique_seq)):
613
614 if sorted_uni_arms[i][0] in unique_seq[j][2].split("_")[0]:
615
616 temp.append(unique_seq[j])
617
618 temp = sorted(temp, key=lambda x: int(x[0].split('-')[1]), reverse=True)
619 fp.write("*********************************************************************************************************\n")
620 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]),"|"))
621 fp.write("*********************************************************************************************************\n\n")
622 [fp.write("%-8s\t%-40s\t%s\n" % (x[0].split("-")[1], x[2],x[9])) for x in temp]
623 fp.write("\n" + "\n")
624 fp.close()
625
626 if f==2:
627
628 if con=="c":
629 fp = open('split3/'+name, 'w')
630 fp.write("%s\t%-42s\t%s\n\n" % ("Number of Reads","Name of isomir","Sequence"))
631 if con=="t":
632 fp = open('split4/'+name, 'w')
633 fp.write("%s\t%-42s\t%s\n\n" % ("Number of Reads","Name of isomir","Sequence"))
634
635
636 for i in range(len(sorted_uni_arms)):
637 temp = []
638 for j in range(len(unique_seq)):
639 if sorted_uni_arms[i][0]==unique_seq[j][2].split("__")[0]:
640 temp.append(unique_seq[j])
641 if temp!=[]:
642 temp = sorted(temp, key=lambda x: int(x[0].split('-')[1]), reverse=True)
643 fp.write("*********************************************************************************************************\n")
644 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]),"|"))
645 fp.write("*********************************************************************************************************\n\n")
646 [fp.write("%-8s\t%-40s\t%s\n" % (x[0].split("-")[1], x[2],x[9])) for x in temp]
647 fp.write("\n" + "\n")
648 fp.close()
649
650
651 ##########################################################################################################################
652
653 def new_mat_seq(pre_unique_seq,mat_mirnas,l):
654
655 unique_iso = []
656 for x in pre_unique_seq:
657 if len(x[2].split("_"))==3:
658 for y in pre_unique_seq:
659 if x[2] in y[2] and int(x[0].split("-")[1])<int(y[0].split("-")[1]):
660 if any(y[2] in lst2 for lst2 in unique_iso)==False:
661 y[2]=">"+y[2]
662 unique_iso.append(y)
663 l.acquire()
664 for x in unique_iso:
665 mat_mirnas.append(x[2])
666 mat_mirnas.append(x[9])
667 l.release()
668
669 #########################################################################################################################
670 def pie_non_temp(merge_LH2E,merge_non_LH2E,merge_LH8E,merge_non_LH8E,c_unmap,t_unmap,c_unmap_counts,t_unmap_counts):
671
672 c_samples=[[x[0],x[1],sum(int(i) for i in x[2:])] for x in merge_LH2E]
673 t_samples=[[x[0],x[1],sum(int(i) for i in x[2:])] for x in merge_LH8E]
674 c_non_samples=[[x[0],x[1],sum(int(i) for i in x[2:])] for x in merge_non_LH2E]
675 t_non_samples=[[x[0],x[1],sum(int(i) for i in x[2:])] for x in merge_non_LH8E]
676
677 c_templ = 0
678 c_tem_counts = 0
679 c_mature = 0
680 c_mat_counts = 0
681 t_templ = 0
682 t_tem_counts = 0
683 t_mature = 0
684 t_mat_counts = 0
685
686 c_non = len(c_non_samples)
687 c_non_counts = sum(x[2] for x in c_non_samples)
688 t_non = len(t_non_samples)
689 t_non_counts = sum(x[2] for x in t_non_samples)
690
691 c_unmap = c_unmap - c_non
692 t_unmap = c_unmap - t_non
693
694 c_unmap_counts=c_unmap_counts - c_non_counts
695 t_unmap_counts=t_unmap_counts - t_non_counts
696
697
698 for x in c_samples:
699
700 if "/" not in x[0]:
701 if "chr" in x[0].split("_")[-1]:
702 c_mature+=1
703 c_mat_counts += x[2]
704 else:
705 c_templ+=1
706 c_tem_counts += x[2]
707 else:
708 f=0
709 for y in x[0].split("/"):
710 if "chr" in y.split("_")[-1]:
711 c_mature+=1
712 c_mat_counts += x[2]
713 f=1
714 break
715 if f==0:
716 c_templ+=1
717 c_tem_counts += x[2]
718
719 for x in t_samples:
720
721 if "/" not in x[0]:
722 if "chr" in x[0].split("_")[-1]:
723 t_mature+=1
724 t_mat_counts += x[2]
725 else:
726 t_templ+=1
727 t_tem_counts += x[2]
728 else:
729 f=0
730 for y in x[0].split("/"):
731 if "chr" in y.split("_")[-1]:
732 t_mature+=1
733 t_mat_counts += x[2]
734 f=1
735 break
736 if f==0:
737 t_templ+=1
738 t_tem_counts += x[2]
739
740 fig = plt.figure(figsize=(7,5))
741 labels = 'miRNA RefSeq','Template', 'Unmapped','Non-template'
742 sizes = [c_mat_counts, c_tem_counts, c_unmap_counts,c_non_counts]
743 colors = ['gold', 'yellowgreen', 'lightcoral', 'lightskyblue']
744 ax1 = plt.subplot2grid((1,2),(0,0))
745 patches, texts, autotexts=plt.pie(sizes, labels=labels, colors=colors, startangle=140,autopct='%1.1f%%',radius=0.8)
746 [x.set_fontsize(8) for x in texts]
747 plt.title('Control Group (reads)',fontsize=12)
748 labels = 'miRNA RefSeq','Template', 'Unmapped','non-template'
749 sizes = [t_mat_counts, t_tem_counts, t_unmap_counts, t_non_counts]
750 colors = ['gold', 'yellowgreen', 'lightcoral', 'lightskyblue']
751 ax2 = plt.subplot2grid((1,2),(0,1))
752 patches, texts, autotexts=plt.pie(sizes, labels=labels, colors=colors, startangle=140,autopct='%1.1f%%',radius=0.8)
753 [x.set_fontsize(8) for x in texts]
754 plt.title('Treated Group (reads)', fontsize=12)
755 plt.savefig('pie_non.png',dpi=300)
756
757 ######################################################################################################################################################
758
759 def merging_names(LH2E_copy,new):
760
761 dupes=[]
762 final_LH2E =[]
763
764 for num in range(len(LH2E_copy)):
765
766 if LH2E_copy[num][1] not in final_LH2E and LH2E_copy[num][0] not in final_LH2E:
767 final_LH2E.append(LH2E_copy[num][1])
768 final_LH2E.append(LH2E_copy[num][0])
769 else:
770 dupes.append(LH2E_copy[num][1])
771
772 dupes=list(set(dupes))
773
774 for i in range(len(dupes)):
775 dupes[i]=[dupes[i]]
776
777 for x in LH2E_copy:
778 for y in dupes:
779 if x[1]==y[0]:
780 fl=0
781 if len(y)==1:
782 y.append(x[0])
783 else:
784 for i in range(1,len(y)):
785 if y[i].split("_")[0]==x[0].split("_")[0]:
786 fl=1
787 if len(x[0])<len(y[i]):
788 del y[i]
789 y.append(x[0])
790 break
791
792 if fl==0:
793 y.append((x[0]))
794
795 for y in dupes:
796 if len(y)>2:
797 for i in range(len(y)-1,1,-1):
798 y[1]=y[1]+"/"+y[i]
799 del y[i]
800
801
802 for x in LH2E_copy:
803 for y in dupes:
804 if x[1]==y[0]:
805 x[0]=y[1]
806
807 LH2E_copy.sort()
808 LH2E_copy=list(LH2E_copy for LH2E_copy,_ in itertools.groupby(LH2E_copy))
809
810 new.extend(LH2E_copy)
811
812
813 ######################################################################################################################################################
814 def pie_temp(merge_LH2E,c_unmap,c_unmap_counts,merge_LH8E,t_unmap,t_unmap_counts):
815
816 c_samples=[[x[0],x[1],sum(int(i) for i in x[2:])] for x in merge_LH2E]
817 t_samples=[[x[0],x[1],sum(int(i) for i in x[2:])] for x in merge_LH8E]
818
819 c_templ = 0
820 c_tem_counts = 0
821 c_mature = 0
822 c_mat_counts = 0
823 t_templ = 0
824 t_tem_counts = 0
825 t_mature = 0
826 t_mat_counts = 0
827
828 for x in c_samples:
829
830 if "/" not in x[0]:
831 if "chr" in x[0].split("_")[-1]:
832 c_mature+=1
833 c_mat_counts += x[2]
834 else:
835 c_templ+=1
836 c_tem_counts += x[2]
837 else:
838 f=0
839 for y in x[0].split("/"):
840 if "chr" in y.split("_")[-1]:
841 c_mature+=1
842 c_mat_counts += x[2]
843 f=1
844 break
845 if f==0:
846 c_templ+=1
847 c_tem_counts += x[2]
848
849 for x in t_samples:
850
851 if "/" not in x[0]:
852 if "chr" in x[0].split("_")[-1]:
853 t_mature+=1
854 t_mat_counts += x[2]
855 else:
856 t_templ+=1
857 t_tem_counts += x[2]
858 else:
859 f=0
860 for y in x[0].split("/"):
861 if "chr" in y.split("_")[-1]:
862 t_mature+=1
863 t_mat_counts += x[2]
864 f=1
865 break
866 if f==0:
867 t_templ+=1
868 t_tem_counts += x[2]
869
870
871 fig = plt.figure()
872 labels = 'miRNA RefSeq','Template', 'Unmapped'
873 sizes = [c_mat_counts, c_tem_counts, c_unmap_counts]
874 colors = ['gold', 'yellowgreen', 'lightskyblue']
875 explode = (0.2, 0.05, 0.1)
876 ax1 = plt.subplot2grid((1,2),(0,0))
877 patches, texts, autotexts=plt.pie(sizes, labels=labels, colors=colors, startangle=140,autopct='%1.1f%%',radius=0.8)
878 [x.set_fontsize(8) for x in texts]
879 plt.title('Control group (reads)', fontsize=12)
880 labels = 'miRNA RefSeq','Template', 'Unmapped'
881 sizes = [t_mat_counts, t_tem_counts, t_unmap_counts]
882 colors = ['gold', 'yellowgreen', 'lightskyblue']
883 explode = (0.2, 0.05, 0.1)
884 ax2 = plt.subplot2grid((1,2),(0,1))
885 patches, texts, autotexts=plt.pie(sizes, labels=labels, colors=colors, startangle=140,autopct='%1.1f%%',radius=0.8)
886 [x.set_fontsize(8) for x in texts]
887 plt.title('Treated group (reads)',fontsize = 12)
888 plt.savefig('pie_tem.png',dpi=300)
889
890 ###################################################################################################################################################################################################################
891
892 def make_spider(merge_LH2E,merge_LH8E):
893
894 c_samples=[[x[0],x[1],sum(int(i) for i in x[2:])] for x in merge_LH2E]
895 t_samples=[[x[0],x[1],sum(int(i) for i in x[2:])] for x in merge_LH8E]
896
897 c_5 = 0
898 c_5_counts = 0
899 c_3 = 0
900 c_3_counts = 0
901 c_both =0
902 c_both_counts=0
903 c_mature = 0
904 c_mat_counts = 0
905 c_exception=0
906 c_exception_counts=0
907
908
909 t_5 = 0
910 t_5_counts = 0
911 t_3 = 0
912 t_3_counts = 0
913 t_both = 0
914 t_both_counts = 0
915 t_mature = 0
916 t_mat_counts = 0
917 t_exception = 0
918 t_exception_counts=0
919
920 for x in c_samples:
921
922 if "/" not in x[0]:
923 if "chr" in x[0].split("_")[-1]:
924 c_mature+=1
925 c_mat_counts += x[2]
926 elif 0 == int(x[0].split("_")[-1]):
927 c_5+=1
928 c_5_counts += x[2]
929 elif 0 == int(x[0].split("_")[-2]):
930 c_3+=1
931 c_3_counts += x[2]
932 else:
933 c_both+=1
934 c_both_counts+=x[2]
935
936 else:
937 f=0
938 for y in x[0].split("/"):
939 if "chr" in y.split("_")[-1]:
940 c_mature+=1
941 c_mat_counts += x[2]
942 f=1
943 break
944 if f==0:
945 for y in x[0].split("/"):
946 c_exception+=1
947 c_exception_counts += x[2]
948
949
950 for x in t_samples:
951
952 if "/" not in x[0]:
953 if "chr" in x[0].split("_")[-1]:
954 t_mature+=1
955 t_mat_counts += x[2]
956 elif 0 == int(x[0].split("_")[-1]):
957 t_5+=1
958 t_5_counts += x[2]
959 elif 0 == int(x[0].split("_")[-2]):
960 t_3+=1
961 t_3_counts += x[2]
962 else:
963 t_both+=1
964 t_both_counts+=x[2]
965
966 else:
967 f=0
968 for y in x[0].split("/"):
969 if "chr" in y.split("_")[-1]:
970 t_mature+=1
971 t_mat_counts += x[2]
972 f=1
973 break
974 if f==0:
975 for y in x[0].split("/"):
976 t_exception+=1
977 t_exception_counts += x[2]
978
979
980 c_all = c_5+c_3+c_both+c_mature+c_exception
981 c_all_counts = c_5_counts + c_3_counts + c_both_counts + c_mat_counts + c_exception_counts
982
983 t_all = t_5+t_3+t_both+t_mature + t_exception
984 t_all_counts = t_5_counts + t_3_counts + t_both_counts + t_mat_counts + t_exception_counts
985
986 c_5 = round(c_5/c_all*100,2)
987 c_3 = round(c_3/c_all*100,2)
988 c_both = round(c_both/c_all*100,2)
989 c_mature = round(c_mature/c_all*100,2)
990 c_exception = round(c_exception/c_all*100,2)
991
992 c_5_counts = round(c_5_counts/c_all_counts*100,2)
993 c_3_counts = round(c_3_counts/c_all_counts*100,2)
994 c_both_counts = round(c_both_counts/c_all_counts*100,2)
995 c_mat_counts = round(c_mat_counts/c_all_counts*100,2)
996 c_exception_counts = round(c_exception_counts/c_all_counts*100,2)
997
998 t_5 = round(t_5/t_all*100,2)
999 t_3 = round(t_3/t_all*100,2)
1000 t_both = round(t_both/t_all*100,2)
1001 t_mature = round(t_mature/t_all*100,2)
1002 t_exception = round(t_exception/t_all*100,2)
1003
1004 t_5_counts = round(t_5_counts/t_all_counts*100,2)
1005 t_3_counts = round(t_3_counts/t_all_counts*100,2)
1006 t_both_counts = round(t_both_counts/t_all_counts*100,2)
1007 t_mat_counts = round(t_mat_counts/t_all_counts*100,2)
1008 t_exception_counts = round(t_exception_counts/t_all_counts*100,2)
1009
1010 radar_max = max(c_5, c_3, c_both,c_mature,c_exception,t_5,t_3,t_both,t_mature,t_exception)
1011 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)
1012
1013 df=pd.DataFrame({
1014 'group':['Controls','Treated'],
1015 """5' and 3' isomiRs""":[c_both,t_both],
1016 """3' isomiRs""":[c_3,t_3],
1017 'miRNA RefSeq':[c_mature,t_mature],
1018 """5' isomiRs""":[c_5,t_5],
1019 'Others*':[c_exception,t_exception]})
1020
1021 df1=pd.DataFrame({
1022 'group':['Controls','Treated'],
1023 """5' and 3' isomiRs""":[c_both_counts,t_both_counts],
1024 """3' isomiRs""":[c_3_counts,t_3_counts],
1025 'miRNA RefSeq':[c_mat_counts,t_mat_counts],
1026 """5' isomiRs""":[c_5_counts,t_5_counts],
1027 'Others*':[c_exception_counts,t_exception_counts]})
1028
1029 spider_last(df,radar_max,1)
1030 spider_last(df1,radar_max_counts,2)
1031
1032
1033 def spider_last(df,radar_max,flag):
1034 # ------- PART 1: Create background
1035 fig = plt.figure()
1036 # number of variable
1037 categories=list(df)[1:]
1038 N = len(categories)
1039
1040 # What will be the angle of each axis in the plot? (we divide the plot / number of variable)
1041 angles = [n / float(N) * 2 * pi for n in range(N)]
1042 angles += angles[:1]
1043
1044 # Initialise the spider plot
1045 ax = plt.subplot(111, polar=True)
1046
1047 # If you want the first axis to be on top:
1048 ax.set_theta_offset(pi/2)
1049 ax.set_theta_direction(-1)
1050
1051 # Draw one axe per variable + add labels labels yet
1052 plt.xticks(angles[:-1], categories, fontsize=11)
1053
1054 # Draw ylabels
1055 radar_max=round(radar_max+radar_max*0.1)
1056 mul=len(str(radar_max))-1
1057 maxi=int(math.ceil(radar_max / pow(10,mul))) * pow(10,mul)
1058 sep = round(maxi/4)
1059 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)
1060 plt.ylim(0, maxi)
1061
1062 # ------- PART 2: Add plots
1063
1064 # Plot each individual = each line of the data
1065 # I don't do a loop, because plotting more than 3 groups makes the chart unreadable
1066
1067 # Ind1
1068 values=df.loc[0].drop('group').values.flatten().tolist()
1069 values += values[:1]
1070 ax.plot(angles, values,'-o', linewidth=1, linestyle='solid', label="Controls")
1071 ax.fill(angles, values, 'b', alpha=0.1)
1072
1073 # Ind2
1074 values=df.loc[1].drop('group').values.flatten().tolist()
1075 values += values[:1]
1076 ax.plot(angles, values, '-o' ,linewidth=1, linestyle='solid', label="Treated")
1077 ax.fill(angles, values, 'r', alpha=0.1)
1078
1079 # Add legend
1080 if flag==1:
1081 plt.legend(loc='upper right', bbox_to_anchor=(0.0, 0.1))
1082 plt.savefig('spider_non_red.png',dpi=300)
1083 else:
1084 plt.legend(loc='upper right', bbox_to_anchor=(0.0, 0.1))
1085 plt.savefig('spider_red.png',dpi=300)
1086
1087
1088 #############################################################################################################################################################################################################
1089 def hist_red(samples,flag):
1090 lengths=[]
1091 cat=[]
1092 total_reads=0
1093 seq=[]
1094
1095 if flag == "c":
1096 title = "Length Distribution of Control group (Redudant reads)"
1097 if flag == "t":
1098 title = "Length Distribution of Treated group (Redudant reads)"
1099
1100 for i in samples:
1101 for x in i:
1102 lengths.append(len(x[9]))
1103 if x[1]=="0":
1104 seq.append([x[9],x[0].split("-")[1],"Mapped"])
1105 cat.append("Mapped")
1106 if x[1] == "4":
1107 seq.append([x[9],x[0].split("-")[1],"Unmapped"])
1108 cat.append("Unmapped")
1109
1110 uni_len=list(set(lengths))
1111 uni_len=[x for x in uni_len if x<=35]
1112 low=min(lengths)
1113 up=max(lengths)
1114 seq.sort()
1115 uni_seq=list(seq for seq,_ in itertools.groupby(seq))
1116 dim=up-low
1117
1118 if dim>20:
1119 s=5
1120 else:
1121 s=8
1122
1123 total_reads+=sum([int(x[1]) for x in uni_seq])
1124
1125 map_reads=[]
1126 unmap_reads=[]
1127 length=[]
1128 for y in uni_len:
1129 map_temp=0
1130 unmap_temp=0
1131 for x in uni_seq:
1132 if len(x[0])==y and x[2]=="Mapped":
1133 map_temp+=int(x[1])
1134 if len(x[0])==y and x[2]=="Unmapped":
1135 unmap_temp+=int(x[1])
1136 if y<=35:
1137 length.append(y)
1138 map_reads.append(round(map_temp/total_reads*100,2))
1139 unmap_reads.append(round(unmap_temp/total_reads*100,2))
1140
1141 ylim=max([sum(x) for x in zip(unmap_reads, map_reads)])
1142 ylim=ylim+ylim*20/100
1143 fig, ax = plt.subplots()
1144 width=0.8
1145 ax.bar(length, unmap_reads, width, label='Unmapped')
1146 h=ax.bar(length, map_reads, width, bottom=unmap_reads, label='Mapped')
1147 plt.xticks(np.arange(length[0], length[-1]+1, 1))
1148 plt.xlabel('Length (nt)',fontsize=14)
1149 plt.ylabel('Percentage',fontsize=14)
1150 plt.title(title,fontsize=14)
1151 ax.legend()
1152 plt.ylim([0, ylim])
1153 ax.grid(axis='y',linewidth=0.2)
1154
1155 if flag=='c':
1156 plt.savefig('c_hist_red.png',dpi=300)
1157
1158 if flag=='t':
1159 plt.savefig('t_hist_red.png',dpi=300)
1160
1161 #################################################################################################################
1162
1163
1164 def logo_seq_red(merge, flag):
1165
1166 if flag=="c":
1167 titlos="Control group (Redundant)"
1168 file_logo="c_logo.png"
1169 file_bar="c_bar.png"
1170 if flag=="t":
1171 titlos="Treated group (Redundant)"
1172 file_logo="t_logo.png"
1173 file_bar="t_bar.png"
1174
1175 c_samples=[[x[0],x[1],sum(int(i) for i in x[2:])] for x in merge]
1176
1177 A=[0]*5
1178 C=[0]*5
1179 G=[0]*5
1180 T=[0]*5
1181 total_reads=0
1182
1183 for y in c_samples:
1184 if "/" in y[0]:
1185 length=[]
1186 for x in y[0].split("/"):
1187 length.append([len(x.split("_")[-1]),x.split("_")[-1],y[2]])
1188
1189 best=min(length)
1190 total_reads+=best[2]
1191 for i in range(5):
1192 if i<len(best[1]):
1193 if best[1][i] == "A":
1194 A[i]+=best[2]
1195 elif best[1][i] == "C":
1196 C[i]+=best[2]
1197 elif best[1][i] == "G":
1198 G[i]+=best[2]
1199 else:
1200 T[i]+=best[2]
1201 else:
1202 total_reads+=y[2]
1203 for i in range(5):
1204 if i<len(y[0].split("_")[-1]):
1205 if y[0].split("_")[-1][i] == "A":
1206 A[i]+=(y[2])
1207 elif y[0].split("_")[-1][i] == "C":
1208 C[i]+=(y[2])
1209 elif y[0].split("_")[-1][i] == "G":
1210 G[i]+=(y[2])
1211 else:
1212 T[i]+=y[2]
1213
1214 A[:] = [round(x*100,1) / total_reads for x in A]
1215 C[:] = [round(x*100,1) / total_reads for x in C]
1216 G[:] = [round(x*100,1) / total_reads for x in G]
1217 T[:] = [round(x*100,1) / total_reads for x in T]
1218
1219
1220
1221 data = {'A':A,'C':C,'G':G,'T':T}
1222 df = pd.DataFrame(data, index=[1,2,3,4,5])
1223 h=df.plot.bar(color=tuple(["g", "b","gold","r"]) )
1224 h.grid(axis='y',linewidth=0.2)
1225 plt.xticks(rotation=0, ha="right")
1226 plt.ylabel("Counts (%)",fontsize=18)
1227 plt.xlabel("Positions (nt)",fontsize=18)
1228 plt.title(titlos,fontsize=20)
1229 plt.tight_layout()
1230 plt.savefig(file_bar, dpi=300)
1231
1232 import logomaker as lm
1233 crp_logo = lm.Logo(df, font_name = 'DejaVu Sans')
1234 crp_logo.style_spines(visible=False)
1235 crp_logo.style_spines(spines=['left', 'bottom'], visible=True)
1236 crp_logo.style_xticks(rotation=0, fmt='%d', anchor=0)
1237
1238 # style using Axes methods
1239 crp_logo.ax.set_title(titlos,fontsize=18)
1240 crp_logo.ax.set_ylabel("Counts (%)", fontsize=16,labelpad=5)
1241 crp_logo.ax.set_xlabel("Positions (nt)",fontsize=16, labelpad=5)
1242 crp_logo.ax.xaxis.set_ticks_position('none')
1243 crp_logo.ax.xaxis.set_tick_params(pad=-1)
1244 figure = plt.gcf()
1245 figure.set_size_inches(6, 4)
1246 crp_logo.fig.savefig(file_logo,dpi=300)
1247
1248 ##########################################################################################################################################################################################################
1249
1250
1251
1252 def logo_seq_non_red(merge_LH2E):
1253
1254 c_samples=[[x[0],x[1],sum(int(i) for i in x[2:])] for x in merge_LH2E]
1255
1256 A=[0]*5
1257 C=[0]*5
1258 G=[0]*5
1259 T=[0]*5
1260
1261 for y in c_samples:
1262 if "/" in y[0]:
1263 length=[]
1264 for x in y[0].split("/"):
1265 length.append([len(x.split("_")[-1]),x.split("_")[-1],y[2]])
1266
1267 best=min(length)
1268 for i in range(5):
1269 if i<len(best[1]):
1270 if best[1][i] == "A":
1271 A[i]+=1
1272 elif best[1][i] == "C":
1273 C[i]+=1
1274 elif best[1][i] == "G":
1275 G[i]+=1
1276 else:
1277 T[i]+=1
1278 else:
1279 for i in range(5):
1280 if i<len(y[0].split("_")[-1]):
1281 if y[0].split("_")[-1][i] == "A":
1282 A[i]+=1
1283 elif y[0].split("_")[-1][i] == "C":
1284 C[i]+=1
1285 elif y[0].split("_")[-1][i] == "G":
1286 G[i]+=1
1287 else:
1288 T[i]+=1
1289
1290 data = {'A':A,'C':C,'G':G,'T':T}
1291 df = pd.DataFrame(data, index=[1,2,3,4,5])
1292 h=df.plot.bar(title="Non-templated nucleotides after templated sequence",color=tuple(["g", "b","gold","r"]))
1293 h.set_xlabel("Positions (nt)")
1294 h.set_ylabel("Unique sequences")
1295 plt.xticks(rotation=0, ha="right")
1296 plt.tight_layout()
1297 plt.savefig("bar2.png", dpi=300)
1298
1299
1300 import logomaker as lm
1301 crp_logo = lm.Logo(df, font_name = 'DejaVu Sans')
1302
1303 # style using Logo methods
1304 crp_logo.style_spines(visible=False)
1305 crp_logo.style_spines(spines=['left', 'bottom'], visible=True)
1306 crp_logo.style_xticks(rotation=0, fmt='%d', anchor=0)
1307
1308 # style using Axes methods
1309 crp_logo.ax.set_ylabel("Unique sequences", labelpad=5)
1310 crp_logo.ax.set_xlabel("Positions (nt)", labelpad=5)
1311 crp_logo.ax.xaxis.set_ticks_position('none')
1312 crp_logo.ax.xaxis.set_tick_params(pad=-1)
1313 crp_logo.ax.set_title("Non-redundant")
1314 figure = plt.gcf()
1315 crp_logo.fig.savefig('logo2.png', dpi=300)
1316
1317
1318 ###################################################################################################################################################################################################################
1319
1320 def ssamples1(tem_names,tem_samp,non_names,non_samp,folder,pro):
1321
1322 for i in range(2,len(tem_samp[0])):
1323
1324 fp = open(folder+tem_names[i-2]+'.txt','w')
1325 fp.write("miRNA id"+"\t"+tem_names[i-2]+"\n")
1326
1327 for x in tem_samp:
1328 fp.write("%s" % "\t".join([x[0],x[i]])+"\n")
1329
1330 for j in range(len(non_names)):
1331 if non_names[j]==tem_names[i-2]:
1332 for x in non_samp:
1333 fp.write("%s" % "\t".join([x[0],x[j+2]])+"\n")
1334 fp.close()
1335
1336 ###################################################################################################################################################################################################################
1337
1338 def download_matures(matures,org_name):
1339
1340 #url = 'ftp://mirbase.org/pub/mirbase/21/mature.fa.gz'
1341 url = 'ftp://mirbase.org/pub/mirbase/CURRENT/mature.fa.gz'
1342 data = urllib.request.urlopen(url).read()
1343 file_mirna = gzip.decompress(data).decode('utf-8')
1344 file_mirna = file_mirna.split("\n")
1345
1346 for i in range(0,len(file_mirna)-1,2):
1347
1348 if org_name in file_mirna[i]:
1349 matures.append(file_mirna[i])
1350 matures.append(file_mirna[i+1])
1351
1352 ###################################################################################################################################################################################################################
1353 def non_template_ref(sc,st,all_isoforms):
1354
1355 pre_uni_seq_con = list(sc)
1356 pre_uni_seq_tre = list(st)
1357
1358 for x in pre_uni_seq_con:
1359 for y in x:
1360 if ">"+y[2] not in all_isoforms and ")_" in y[2] :
1361 all_isoforms.append(">"+y[2])
1362 all_isoforms.append(y[9])
1363
1364
1365 for x in pre_uni_seq_tre:
1366 for y in x:
1367 if ">"+y[2] not in all_isoforms and ")_" in y[2]:
1368 all_isoforms.append(">"+y[2])
1369 all_isoforms.append(y[9])
1370
1371 ################################################################################################################################################################################################
1372
1373 def deseqe2(sample,mir_names,l,new_d,sample_name,sample_order):
1374
1375 for y in mir_names:
1376 flag=0
1377 for x in sample:
1378 if y[0]==x[0]:
1379 flag=1
1380 break
1381 if flag==0:
1382 sample.append([y[0],"0",y[1]])
1383
1384 sample.sort(key=lambda x: x[0])
1385 sample=list(sample for sample,_ in itertools.groupby(sample))
1386
1387 l.acquire()
1388 new_d.append(sample)
1389 sample_order.append(sample_name)
1390 l.release()
1391
1392 ###############################################################################################################################################################################################
1393
1394 if __name__ == '__main__':
1395
1396 starttime = time.time()
1397
1398 q1 = Queue()
1399 q2 = Queue()
1400 lock = Lock()
1401 manager = Manager()
1402
1403 mature_mirnas=manager.list()
1404 ps_mature=Process(target=download_matures,args=(mature_mirnas,args.org_name))
1405 ps_mature.start()
1406
1407 args.control[0]=args.control[0][1:]
1408 args.control[len(args.control)-1][:-1]
1409 control = [(args.control[i:i+2]) for i in range(0, len(args.control), 2)]
1410
1411 args.treated[0]=args.treated[0][1:]
1412 args.treated[len(args.treated)-1][:-1]
1413 treated = [(args.treated[i:i+2]) for i in range(0, len(args.treated), 2)]
1414
1415
1416 ############## Detection of templated isoforms ################
1417
1418 radar = manager.list([0,0,0,0])
1419 samples = manager.list()
1420 data= manager.list()
1421 names_con=manager.list()
1422 samples_mirna_names=manager.list()
1423 deseq=manager.list()
1424 unmap_seq=manager.Value('i',0)
1425 unmap_counts=manager.Value('i',0)
1426 LH2E_names=manager.list()
1427 ini_c_samples = manager.list()
1428
1429
1430 radar1 = manager.list([0,0,0,0])
1431 samples1 = manager.list()
1432 data1 = manager.list()
1433 names_tre = manager.list()
1434 samples_mirna_names1=manager.list()
1435 deseq1=manager.list()
1436 unmap1_seq = manager.Value('i',0)
1437 unmap1_counts = manager.Value('i',0)
1438 LH8E_names=manager.list()
1439 ini_t_samples = manager.list()
1440 ps_mature.join()
1441
1442
1443 mature_mirnas=list(mature_mirnas)
1444
1445
1446 starttime1 = time.time()
1447 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]
1448 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])
1449
1450 [p.start() for p in ps_sam]
1451 [p.join() for p in ps_sam]
1452 print('SAM took {} seconds'.format(time.time() - starttime1))
1453
1454 ps_hist=[Process(target=hist_red,args=(ini_c_samples,'c'))]
1455 ps_hist.extend([Process(target=hist_red,args=(ini_t_samples,'t'))])
1456 [x.start() for x in ps_hist]
1457
1458 starttime200=time.time()
1459
1460 sc = list(samples)
1461 st = list(samples1)
1462
1463 names_con=list(names_con)
1464 names_tre=list(names_tre)
1465 samples_mirna_names=list(samples_mirna_names)
1466 samples_mirna_names.sort()
1467 samples_mirna_names=list(samples_mirna_names for samples_mirna_names,_ in itertools.groupby(samples_mirna_names))
1468
1469 samples_mirna_names1=list(samples_mirna_names1)
1470 samples_mirna_names1.sort()
1471 samples_mirna_names1=list(samples_mirna_names1 for samples_mirna_names1,_ in itertools.groupby(samples_mirna_names1))
1472
1473 deseq=list(deseq)
1474 deseq1=list(deseq1)
1475
1476 new_names_con=manager.list()
1477 new_names_tre=manager.list()
1478 new_deseq=manager.list()
1479 new_deseq1=manager.list()
1480 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)]
1481 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)])
1482
1483 [z.start() for z in ps_deseq]
1484 [z.join() for z in ps_deseq]
1485 new_deseq=list(new_deseq)
1486 new_deseq1=list(new_deseq1)
1487
1488 LH2E=[[x[0],x[2]] for x in new_deseq[0]]
1489 [LH2E[i].append(y[i][1]) for i,_ in enumerate(LH2E) for y in new_deseq]
1490
1491 LH8E=[[x[0],x[2]] for x in new_deseq1[0]]
1492 [LH8E[i].append(y[i][1]) for i,_ in enumerate(LH8E) for y in new_deseq1]
1493
1494 print('Deseq took {} seconds'.format(time.time() - starttime200))
1495
1496 merg_nam_LH2E=manager.list()
1497 merg_nam_LH8E=manager.list()
1498
1499 LH2E_copy=copy.deepcopy(list(LH2E))
1500 LH8E_copy=copy.deepcopy(list(LH8E))
1501
1502 fil_sort_tre=manager.list()
1503 fil_sort_con=manager.list()
1504 raw_sort_tre=manager.list()
1505 raw_sort_con=manager.list()
1506
1507 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))
1508 ps_main.start()
1509
1510 if args.anal=="2":
1511 all_iso = manager.list()
1512 ps_non_iso = Process(target=non_template_ref,args=(sc,st,all_iso))
1513 ps_non_iso.start()
1514
1515 ps_merge = [Process(target=merging_names,args=(LH2E_copy,merg_nam_LH2E))]
1516 ps_merge.extend([Process(target=merging_names,args=(LH8E_copy,merg_nam_LH8E))])
1517 [x.start() for x in ps_merge]
1518 [x.join() for x in ps_merge]
1519
1520 merg_nam_LH2E=list(merg_nam_LH2E)
1521 merg_nam_LH8E=list(merg_nam_LH8E)
1522
1523 starttime2 = time.time()
1524 procs = [Process(target=DB_write,args=(x[0],x[1],x[2],x[3],1)) for x in data]
1525 procs.extend([Process(target=DB_write,args=(x[0],x[1],x[2],x[3],1)) for x in data1])
1526 procs.extend([Process(target=make_spider,args=(merg_nam_LH2E,merg_nam_LH8E))])
1527 if args.anal == "1":
1528 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))])
1529
1530 [p.start() for p in procs]
1531
1532
1533 if args.anal=="1":
1534 [x.join() for x in ps_hist]
1535 [p.join() for p in procs]
1536 ps_pdf = Process(target=pdf_before_DE,args=(args.anal))
1537 ps_pdf.start()
1538
1539 print('Graphs took {} seconds'.format(time.time() - starttime2))
1540
1541 ps_main.join()
1542
1543 fil_sort_con=list(fil_sort_con)
1544 fil_sort_tre=list(fil_sort_tre)
1545 if fil_sort_con==[]:
1546 fil_sort_con=raw_sort_con
1547 fil_sort_tre=raw_sort_tre
1548
1549 raw_sort_con=list(raw_sort_con)
1550 raw_sort_tre=list(raw_sort_tre)
1551 names_con=list(new_names_con)
1552 names_tre=list(new_names_tre)
1553
1554 ps_write = Process(target=write_main,args=(raw_sort_con, raw_sort_tre, fil_sort_con, fil_sort_tre, names_con,names_tre,1))
1555 ps_write.start()
1556
1557 ps1_matrix = [Process(target=ssamples,args=(names_con,fil_sort_con,"Diff/temp_con/",0))]
1558 ps1_matrix.extend([Process(target=ssamples,args=(names_tre,fil_sort_tre,"Diff/temp_tre/",0))])
1559 [p.start() for p in ps1_matrix]
1560
1561 if args.anal=="1":
1562 ps_pdf.join()
1563 if args.anal=="2":
1564 [p.join() for p in procs]
1565 [x.join() for x in ps_hist]
1566
1567 ps_write.join()
1568 [p.join() for p in ps1_matrix]
1569
1570
1571 ############################## Detection of Both #######################################
1572
1573 starttime10 = time.time()
1574
1575 if args.anal == "2":
1576
1577 n_data= manager.list()
1578 n_names_con=manager.list()
1579 n_samples_mirna_names=manager.list()
1580 n_deseq=manager.list()
1581 n_LH2E_names=manager.list()
1582
1583 n_data1 = manager.list()
1584 n_names_tre = manager.list()
1585 n_samples_mirna_names1=manager.list()
1586 n_deseq1=manager.list()
1587 n_LH8E_names=manager.list()
1588
1589
1590 new_mat_mirnas = list(mature_mirnas)
1591 ps_non_iso.join()
1592
1593 all_iso=list(all_iso)
1594 new_mat_mirnas.extend(all_iso)
1595
1596 starttime11=time.time()
1597
1598 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]
1599 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])
1600
1601 [p.start() for p in ps_sam]
1602 [p.join() for p in ps_sam]
1603
1604 print('Non-sam took {} seconds'.format(time.time() - starttime11))
1605
1606 starttime12=time.time()
1607
1608 n_names_con=list(n_names_con)
1609 n_names_tre=list(n_names_tre)
1610 n_samples_mirna_names=list(n_samples_mirna_names)
1611 n_samples_mirna_names.sort()
1612 n_samples_mirna_names=list(n_samples_mirna_names for n_samples_mirna_names,_ in itertools.groupby(n_samples_mirna_names))
1613
1614 n_samples_mirna_names1=list(n_samples_mirna_names1)
1615 n_samples_mirna_names1.sort()
1616 n_samples_mirna_names1=list(n_samples_mirna_names1 for n_samples_mirna_names1,_ in itertools.groupby(n_samples_mirna_names1))
1617
1618 n_deseq=list(n_deseq)
1619 n_deseq1=list(n_deseq1)
1620
1621 new_n_names_con=manager.list()
1622 new_n_names_tre=manager.list()
1623 n_new_deseq=manager.list()
1624 n_new_deseq1=manager.list()
1625 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)]
1626 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)])
1627
1628 [x.start() for x in ps_deseq]
1629 [x.join() for x in ps_deseq]
1630 n_new_deseq=list(n_new_deseq)
1631 n_new_deseq1=list(n_new_deseq1)
1632
1633 n_LH2E=[[x[0],x[2]] for x in n_new_deseq[0]]
1634 [n_LH2E[i].append(y[i][1]) for i,_ in enumerate(n_LH2E) for y in n_new_deseq]
1635
1636 n_LH8E=[[x[0],x[2]] for x in n_new_deseq1[0]]
1637 [n_LH8E[i].append(y[i][1]) for i,_ in enumerate(n_LH8E) for y in n_new_deseq1]
1638
1639 print('Non-deseq took {} seconds'.format(time.time() - starttime12))
1640
1641 merg_nam_n_LH2E=manager.list()
1642 merg_nam_n_LH8E=manager.list()
1643
1644 n_LH2E_copy=copy.deepcopy(list(n_LH2E))
1645 n_LH8E_copy=copy.deepcopy(list(n_LH8E))
1646
1647 n_sort_tre=manager.list()
1648 n_sort_con=manager.list()
1649
1650 n_fil_sort_con=manager.list()
1651 n_fil_sort_tre=manager.list()
1652 n_raw_sort_con=manager.list()
1653 n_raw_sort_tre=manager.list()
1654
1655 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))
1656 ps_main.start()
1657
1658 ps_merge = [Process(target=merging_names,args=(n_LH2E_copy,merg_nam_n_LH2E))]
1659 ps_merge.extend([Process(target=merging_names,args=(n_LH8E_copy,merg_nam_n_LH8E))])
1660 [p.start() for p in ps_merge]
1661 [p.join() for p in ps_merge]
1662
1663 merg_nam_n_LH2E=list(merg_nam_n_LH2E)
1664 merg_nam_n_LH8E=list(merg_nam_n_LH8E)
1665
1666 procs = [Process(target=DB_write,args=(x[0],x[1],x[2],x[3],2)) for x in n_data]
1667 procs.extend([Process(target=DB_write,args=(x[0],x[1],x[2],x[3],2)) for x in n_data1])
1668 procs.extend([Process(target=logo_seq_red,args=(merg_nam_n_LH2E,'c'))])
1669 procs.extend([Process(target=logo_seq_red,args=(merg_nam_n_LH8E,'t'))])
1670 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))])
1671
1672 starttime13=time.time()
1673 [p.start() for p in procs]
1674 [p.join() for p in procs]
1675
1676 print('Graphs took {} seconds'.format(time.time() - starttime13))
1677
1678 procs1 = Process(target=pdf_before_DE,args=(args.anal))
1679 procs1.start()
1680
1681 starttime14=time.time()
1682 ps_main.join()
1683
1684 n_fil_sort_con=list(n_fil_sort_con)
1685 n_fil_sort_tre=list(n_fil_sort_tre)
1686 if n_fil_sort_con==[]:
1687 n_fil_sort_con=n_raw_sort_con
1688 n_fil_sort_tre=n_raw_sort_tre
1689
1690 n_raw_sort_con=list(n_raw_sort_con)
1691 n_raw_sort_tre=list(n_raw_sort_tre)
1692 n_names_con=list(new_n_names_con)
1693 n_names_tre=list(new_n_names_tre)
1694
1695 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))
1696 ps_write.start()
1697
1698 ps1_matrix = [Process(target=ssamples1,args=(n_names_con,n_fil_sort_con,names_con,fil_sort_con,"Diff/n_temp_con/",0))]
1699 ps1_matrix.extend([Process(target=ssamples1,args=(n_names_tre,n_fil_sort_tre,names_tre,fil_sort_tre,"Diff/n_temp_tre/",0))])
1700 [p.start() for p in ps1_matrix]
1701
1702 ps_write.join()
1703 [p.join() for p in ps1_matrix]
1704 procs1.join()
1705 print('That took {} seconds'.format(time.time() - starttime10))
1706 print('That took {} seconds'.format(time.time() - starttime))
1707
1708
1709
1710
1711
1712
1713
1714