|
3
|
1 import itertools
|
|
|
2 import time
|
|
|
3 import sys
|
|
|
4 import os
|
|
|
5 import urllib.request
|
|
|
6 import gzip
|
|
|
7 from multiprocessing import Process, Queue, Lock, Pool, Manager, Value
|
|
|
8 import subprocess
|
|
|
9 import argparse
|
|
|
10 from collections import OrderedDict
|
|
|
11 from matplotlib.backends.backend_pdf import PdfPages
|
|
|
12 import pandas as pd
|
|
|
13 from math import pi
|
|
|
14 import numpy as np
|
|
|
15 import matplotlib.pyplot as plt
|
|
|
16 from matplotlib.ticker import PercentFormatter
|
|
|
17 import seaborn as sns
|
|
|
18 import scipy.stats as stats
|
|
|
19 from plotnine import *
|
|
|
20 import math
|
|
|
21 import re
|
|
|
22 import matplotlib.ticker as mtick
|
|
|
23 import copy
|
|
|
24
|
|
|
25
|
|
|
26
|
|
|
27
|
|
|
28 """---------------------- Simple Functions -----------------------"""
|
|
|
29
|
|
|
30 # Read a file and return it as a list
|
|
|
31 def read(path, flag):
|
|
|
32 if flag == 0:
|
|
|
33 with open(path) as fp:
|
|
|
34 file=fp.readlines()
|
|
|
35 fp.close()
|
|
|
36 return file
|
|
|
37
|
|
|
38 if flag == 1:
|
|
|
39 with open(path) as fp:
|
|
|
40 file = fp.read().splitlines()
|
|
|
41 fp.close()
|
|
|
42 return file
|
|
|
43
|
|
|
44 # Write a list to a txt file
|
|
|
45 def write(path, list):
|
|
|
46 with open(path,'w') as fp:
|
|
|
47 for x in list:
|
|
|
48 fp.write(str("\t".join(x[1:-1])))
|
|
|
49 fp.close()
|
|
|
50
|
|
|
51 """---------------------- RNA-seq Functions ----------------------"""
|
|
|
52
|
|
|
53 # Detect the longest common substring sequence between two mirnas
|
|
|
54 def longestSubstring(str1, str2):
|
|
|
55
|
|
|
56 from difflib import SequenceMatcher
|
|
|
57 # initialize SequenceMatcher object with
|
|
|
58 # input string
|
|
|
59 seqMatch = SequenceMatcher(None, str1, str2)
|
|
|
60
|
|
|
61 # find match of longest sub-string
|
|
|
62 # output will be like Match(a=0, b=0, size=5)
|
|
|
63 match = seqMatch.find_longest_match(0, len(str1), 0, len(str2))
|
|
|
64
|
|
|
65 # print longest substring
|
|
|
66 if (match.size != 0):
|
|
|
67 return str1[match.a: match.a + match.size]
|
|
|
68 else:
|
|
|
69 print('No longest common sub-string found')
|
|
|
70
|
|
|
71
|
|
|
72
|
|
|
73 ########################################################################################################################################################
|
|
|
74 def collapse_sam(path):
|
|
|
75
|
|
|
76 ini_sam=read(path,0)
|
|
|
77 main_sam = [x.rstrip("\n").split("\t") for x in ini_sam if "@" not in x.split("\t")[0]]
|
|
|
78 intro_sam = [x.rstrip("\n").split("\t") for x in ini_sam if "@" in x.split("\t")[0]]
|
|
|
79
|
|
|
80 uni_seq = []
|
|
|
81 for x in main_sam:
|
|
|
82
|
|
|
83 if [x[2], x[9]] not in uni_seq:
|
|
|
84 uni_seq.append([x[2], x[9]])
|
|
|
85
|
|
|
86 new_main_sam=[]
|
|
|
87 incr_num=0
|
|
|
88 for i in range(len(uni_seq)):
|
|
|
89 count=0
|
|
|
90 incr_num+=1
|
|
|
91 for y in main_sam:
|
|
|
92 if uni_seq[i][1]==y[9] and uni_seq[i][0]==y[2]:
|
|
|
93 count+=1
|
|
|
94 temp=y
|
|
|
95 temp[10]="~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
|
|
|
96 temp[0]=str(incr_num)+"-"+str(count)
|
|
|
97 new_main_sam.append(temp)
|
|
|
98
|
|
|
99 new_sam=intro_sam+new_main_sam
|
|
|
100
|
|
|
101 return new_sam
|
|
|
102
|
|
|
103 #################################################################################################################################################################
|
|
|
104
|
|
|
105 def sam(mature_mirnas,path,name,con,l,samples,data,names,unmap_seq,samples_mirna_names,deseq,LHE_names,umi,ini_sample,unmap_counts):
|
|
|
106
|
|
|
107 # read the sam file
|
|
|
108
|
|
|
109 ini_sam=read(path,0)
|
|
|
110 new_main_sam = [x.rstrip("\n").split("\t") for x in ini_sam if "@" not in x.split("\t")[0]]
|
|
|
111 unique_seq = [x for x in new_main_sam if x[1] == '0' and len(x[9])>=18 and len(x[9])<=26]
|
|
|
112
|
|
|
113 # Calculate the shifted positions for every isomir and add them to the name of it
|
|
|
114 sorted_uni_arms = []
|
|
|
115
|
|
|
116 for i in range(len(mature_mirnas)):
|
|
|
117 tmp_count_reads = 0 # calculate the total number of reads
|
|
|
118 tmp_count_seq = 0 # calculate the total number of sequences
|
|
|
119 for j in range(len(unique_seq)):
|
|
|
120
|
|
|
121 if mature_mirnas[i] == unique_seq[j][2]:
|
|
|
122
|
|
|
123 temp_mature = mature_mirnas[i+1]
|
|
|
124 off_part = longestSubstring(temp_mature, unique_seq[j][9])
|
|
|
125
|
|
|
126 mat_diff = temp_mature.split(off_part)
|
|
|
127 mat_diff = [len(mat_diff[0]), len(mat_diff[1])]
|
|
|
128
|
|
|
129 unique_diff = unique_seq[j][9].split(off_part)
|
|
|
130 unique_diff = [len(unique_diff[0]), len(unique_diff[1])]
|
|
|
131
|
|
|
132 # Problem with hsa-miR-8485
|
|
|
133 if mat_diff[1]!=0 and unique_diff[1]!=0:
|
|
|
134 unique_seq[j]=1
|
|
|
135 pre_pos = 0
|
|
|
136 post_pos = 0
|
|
|
137
|
|
|
138 elif mat_diff[0]!=0 and unique_diff[0]!=0:
|
|
|
139 unique_seq[j]=1
|
|
|
140 pre_pos = 0
|
|
|
141 post_pos = 0
|
|
|
142
|
|
|
143 else:
|
|
|
144 pre_pos = mat_diff[0]-unique_diff[0]
|
|
|
145 post_pos = unique_diff[1]-mat_diff[1]
|
|
|
146 tmp_count_reads = tmp_count_reads + int(unique_seq[j][0].split("-")[1])
|
|
|
147 tmp_count_seq = tmp_count_seq+1
|
|
|
148
|
|
|
149
|
|
|
150 if pre_pos != 0 or post_pos != 0:
|
|
|
151 if pre_pos == 0:
|
|
|
152 unique_seq[j][2] = unique_seq[j][2] + "_" +str(pre_pos) + "_" + '{:+d}'.format(post_pos)
|
|
|
153 elif post_pos == 0:
|
|
|
154 unique_seq[j][2] = unique_seq[j][2] + "_" + '{:+d}'.format(pre_pos) + "_" + str(post_pos)
|
|
|
155 else:
|
|
|
156 unique_seq[j][2] = unique_seq[j][2]+"_"+'{:+d}'.format(pre_pos)+"_"+'{:+d}'.format(post_pos)
|
|
|
157
|
|
|
158 for x in range(unique_seq.count(1)):
|
|
|
159 unique_seq.remove(1)
|
|
|
160 if tmp_count_reads != 0 and tmp_count_seq != 0:
|
|
|
161 sorted_uni_arms.append([mature_mirnas[i], tmp_count_seq, tmp_count_reads])
|
|
|
162
|
|
|
163 # Store name of arms, number of sequences and number of reads
|
|
|
164 sorted_uni_arms = sorted(sorted_uni_arms, key=lambda x: x[1], reverse=True)
|
|
|
165
|
|
|
166 for y in sorted_uni_arms:
|
|
|
167 counts=0
|
|
|
168 seqs=0
|
|
|
169 for x in unique_seq:
|
|
|
170 if y[0]==x[2].split("_")[0]+"_"+x[2].split("_")[1]:
|
|
|
171 counts+=int(x[0].split("-")[1])
|
|
|
172 seqs+=1
|
|
|
173
|
|
|
174 y[1]=seqs
|
|
|
175 y[2]=counts
|
|
|
176
|
|
|
177 LHE=[]
|
|
|
178
|
|
|
179 l.acquire()
|
|
|
180 if con=="c":
|
|
|
181 LHE.extend(z[2] for z in unique_seq)
|
|
|
182 for y in unique_seq:
|
|
|
183 samples_mirna_names.append([y[2],y[9]])
|
|
|
184 deseq.append([[x[2], x[0].split('-')[1], x[9]] for x in unique_seq])
|
|
|
185 LHE_names.extend(LHE)
|
|
|
186 unmap_seq.value += sum([1 for x in new_main_sam if x[1] == '4'])
|
|
|
187 unmap_counts.value += sum([int(x[0].split("-")[1]) for x in new_main_sam if x[1] == '4'])
|
|
|
188 names.append(name)
|
|
|
189 samples.append(unique_seq)
|
|
|
190 data.append([con,name,unique_seq,sorted_uni_arms])
|
|
|
191 ini_sample.append(new_main_sam)
|
|
|
192
|
|
|
193 if con=="t":
|
|
|
194 LHE.extend(z[2] for z in unique_seq)
|
|
|
195 for y in unique_seq:
|
|
|
196 samples_mirna_names.append([y[2],y[9]])
|
|
|
197 deseq.append([[x[2], x[0].split('-')[1], x[9]] for x in unique_seq])
|
|
|
198 LHE_names.extend(LHE)
|
|
|
199 unmap_seq.value += sum([1 for x in new_main_sam if x[1] == '4'])
|
|
|
200 unmap_counts.value += sum([int(x[0].split("-")[1]) for x in new_main_sam if x[1] == '4'])
|
|
|
201 names.append(name)
|
|
|
202 samples.append(unique_seq)
|
|
|
203 data.append([con,name,unique_seq,sorted_uni_arms])
|
|
|
204 ini_sample.append(new_main_sam)
|
|
|
205 l.release()
|
|
|
206
|
|
|
207
|
|
|
208 ######################################################################################################################################
|
|
|
209 """
|
|
|
210 Read a sam file from Bowtie and do the followings:
|
|
|
211
|
|
|
212 1) Remove reverse stranded mapped reads
|
|
|
213 2) Remove unmapped reads
|
|
|
214 3) Remove all sequences with reads less than 11 reads
|
|
|
215 4) Sort the arms with the most sequences in decreading rate
|
|
|
216 5) Sort the sequences of every arm with the most reads in decreasing rate
|
|
|
217 6) Calculate total number of sequences of every arm
|
|
|
218 7) Calculate total number of reads of sequences of every arm.
|
|
|
219 8) Store all the informations in a txt file
|
|
|
220
|
|
|
221 """
|
|
|
222
|
|
|
223 def non_sam(mature_mirnas,path,name,con,l,data,names,n_deseq,n_samples_mirna_names,n_LHE_names):
|
|
|
224
|
|
|
225
|
|
|
226 ini_sam=read(path,0)
|
|
|
227 new_main_sam = [x.rstrip("\n").split("\t") for x in ini_sam if "@" not in x.split("\t")[0]]
|
|
|
228 unique_seq=[]
|
|
|
229 unique_seq = [x for x in new_main_sam if x[1] == '4' and len(x[9])>=18 and len(x[9])<=26]
|
|
|
230 uni_seq=[]
|
|
|
231
|
|
|
232 # Calculate the shifted positions for every isomir and add them to the name of it
|
|
|
233 sorted_uni_arms = []
|
|
|
234 for i in range(1,len(mature_mirnas),2):
|
|
|
235 tmp_count_reads = 0 # calculate the total number of reads
|
|
|
236 tmp_count_seq = 0 # calculate the total number of sequences
|
|
|
237
|
|
|
238 for j in range(len(unique_seq)):
|
|
|
239
|
|
|
240 temp_mature = mature_mirnas[i].strip().replace("U", "T")
|
|
|
241
|
|
|
242 if temp_mature in unique_seq[j][9]:
|
|
|
243
|
|
|
244 off_part = longestSubstring(temp_mature, unique_seq[j][9])
|
|
|
245
|
|
|
246 mat_diff = temp_mature.split(off_part)
|
|
|
247 mat_diff = [len(mat_diff[0]), len(mat_diff[1])]
|
|
|
248
|
|
|
249 unique_diff = unique_seq[j][9].split(off_part)
|
|
|
250 if len(unique_diff)<=2:
|
|
|
251 unique_diff = [len(unique_diff[0]), len(unique_diff[1])]
|
|
|
252
|
|
|
253 pre_pos = mat_diff[0]-unique_diff[0]
|
|
|
254 post_pos = unique_diff[1]-mat_diff[1]
|
|
|
255
|
|
|
256 lengthofmir = len(off_part) + post_pos
|
|
|
257 if pre_pos == 0:
|
|
|
258 tmp_count_reads = tmp_count_reads + int(unique_seq[j][0].split("-")[1])
|
|
|
259 tmp_count_seq = tmp_count_seq + 1
|
|
|
260
|
|
|
261 if pre_pos == 0:
|
|
|
262 t_name=copy.deepcopy(unique_seq[j])
|
|
|
263 t_name[2]=mature_mirnas[i - 1] + "__" + str(pre_pos) + "_" + '{:+d}'.format(post_pos) + "_" + str(unique_seq[j][9][len(off_part):])
|
|
|
264 uni_seq.append(t_name)
|
|
|
265
|
|
|
266 if tmp_count_reads != 0 and tmp_count_seq != 0:
|
|
|
267 sorted_uni_arms.append([mature_mirnas[i-1], tmp_count_seq, tmp_count_reads])
|
|
|
268
|
|
|
269
|
|
|
270 sorted_uni_arms = sorted(sorted_uni_arms, key=lambda x: x[1], reverse=True)
|
|
|
271 unique_seq = list(map(list, OrderedDict.fromkeys(map(tuple,uni_seq))))
|
|
|
272
|
|
|
273 LHE=[]
|
|
|
274
|
|
|
275 l.acquire()
|
|
|
276 if con=="c":
|
|
|
277 LHE.extend(x[2] for x in unique_seq if x[2]!="*")
|
|
|
278 for x in unique_seq:
|
|
|
279 if x[2]!="*":
|
|
|
280 n_samples_mirna_names.append([x[2],x[9]])
|
|
|
281 n_deseq.append([[x[2], x[0].split('-')[1], x[9]] for x in unique_seq if x[2]!="*"])
|
|
|
282 n_LHE_names.extend(LHE)
|
|
|
283 names.append(name)
|
|
|
284 data.append([con,name,unique_seq,sorted_uni_arms])
|
|
|
285
|
|
|
286
|
|
|
287 if con=="t":
|
|
|
288 LHE.extend(x[2] for x in unique_seq if x[2]!="*")
|
|
|
289 for x in unique_seq:
|
|
|
290 if x[2]!="*":
|
|
|
291 n_samples_mirna_names.append([x[2],x[9]])
|
|
|
292 n_deseq.append([[x[2], x[0].split('-')[1], x[9]] for x in unique_seq if x[2]!="*"])
|
|
|
293 n_LHE_names.extend(LHE)
|
|
|
294 names.append(name)
|
|
|
295 data.append([con,name,unique_seq,sorted_uni_arms])
|
|
|
296 l.release()
|
|
|
297
|
|
|
298 #################################################################################################################################################################################################################
|
|
|
299
|
|
|
300 def deseq2_temp(samples_mirna_names,deseq,con,l):
|
|
|
301
|
|
|
302 samples_mirna_names.sort(key=lambda x:[0])
|
|
|
303
|
|
|
304 for i in range(len(deseq)):
|
|
|
305 for y in samples_mirna_names:
|
|
|
306 flag = 0
|
|
|
307 for x in deseq[i]:
|
|
|
308 if y[0] == x[0]:
|
|
|
309 flag = 1
|
|
|
310 break
|
|
|
311
|
|
|
312 if flag == 0:
|
|
|
313 deseq[i].append([y[0], "0", y[1]])
|
|
|
314
|
|
|
315 [deseq[i].sort(key=lambda x: x[0]) for i, _ in enumerate(deseq)]
|
|
|
316 deseq_final = [[x[0],x[2]] for x in deseq[0]]
|
|
|
317 [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]]
|
|
|
318
|
|
|
319 l.acquire()
|
|
|
320 if con=="c":
|
|
|
321 q1.put(deseq_final)
|
|
|
322
|
|
|
323 if con=="t":
|
|
|
324 q2.put(deseq_final)
|
|
|
325 l.release()
|
|
|
326
|
|
|
327
|
|
|
328 ##################################################################################################################################################################################################################
|
|
|
329
|
|
|
330 def main_temp(LH2E, LH2E_names, LH8E, LH8E_names,flag,names_con,names_tre,filter_LH8E,filter_LH2E,raw_LH8E,raw_LH2E,per,count):
|
|
|
331
|
|
|
332 LH8E_add_names = [x for x in LH2E_names if x not in LH8E_names]
|
|
|
333 LH2E_add_names = [x for x in LH8E_names if x not in LH2E_names]
|
|
|
334
|
|
|
335 LH8E_add_names.sort()
|
|
|
336 LH2E_add_names.sort()
|
|
|
337 LH8E_add_names = list(LH8E_add_names for LH8E_add_names,_ in itertools.groupby(LH8E_add_names))
|
|
|
338 LH2E_add_names = list(LH2E_add_names for LH2E_add_names,_ in itertools.groupby(LH2E_add_names))
|
|
|
339
|
|
|
340 LH2E.sort()
|
|
|
341 LH8E.sort()
|
|
|
342 LH2E = list(LH2E for LH2E,_ in itertools.groupby(LH2E))
|
|
|
343 LH8E = list(LH8E for LH8E,_ in itertools.groupby(LH8E))
|
|
|
344
|
|
|
345 print("LHE_names")
|
|
|
346 print([len(LH8E_add_names),len(LH2E_add_names)])
|
|
|
347 print([len(LH8E),len(LH2E)])
|
|
|
348
|
|
|
349 zeros=["0"]*(len(LH8E[0])-2)
|
|
|
350 [LH8E_add_names[i].extend(zeros) for i,_ in enumerate(LH8E_add_names)]
|
|
|
351 LH8E=LH8E+LH8E_add_names
|
|
|
352
|
|
|
353 zeros=["0"]*(len(LH2E[0])-2)
|
|
|
354 [LH2E_add_names[i].extend(zeros) for i,_ in enumerate(LH2E_add_names)]
|
|
|
355 LH2E=LH2E+LH2E_add_names
|
|
|
356
|
|
|
357 dupes=[]
|
|
|
358 final_LH2E =[]
|
|
|
359
|
|
|
360 for num,_ in enumerate(LH2E):
|
|
|
361
|
|
|
362 if LH2E[num][1] not in final_LH2E and LH2E[num][0] not in final_LH2E:
|
|
|
363 final_LH2E.append(LH2E[num][1])
|
|
|
364 final_LH2E.append(LH2E[num][0])
|
|
|
365 else:
|
|
|
366 dupes.append(LH2E[num][1])
|
|
|
367
|
|
|
368 dupes=list(set(dupes))
|
|
|
369
|
|
|
370 dupes=[[x] for x in dupes]
|
|
|
371
|
|
|
372 for x in LH2E:
|
|
|
373 for y in dupes:
|
|
|
374 if x[1]==y[0]:
|
|
|
375 fl=0
|
|
|
376 if len(y)==1:
|
|
|
377 y.append(x[0])
|
|
|
378 else:
|
|
|
379 for i in range(1,len(y)):
|
|
|
380 if y[i].split("_")[0]+"_"+y[i].split("_")[1]==x[0].split("_")[0]+"_"+x[0].split("_")[1]:
|
|
|
381 fl=1
|
|
|
382 if len(x[0])<len(y[i]):
|
|
|
383 del y[i]
|
|
|
384 y.append(x[0])
|
|
|
385 break
|
|
|
386
|
|
|
387 if fl==0:
|
|
|
388 y.append((x[0]))
|
|
|
389
|
|
|
390 for y in dupes:
|
|
|
391 if len(y)>2:
|
|
|
392 for i in range(len(y)-1,1,-1):
|
|
|
393 y[1]=y[1]+"/"+y[i]
|
|
|
394 del y[i]
|
|
|
395
|
|
|
396
|
|
|
397 for x in LH2E:
|
|
|
398 for y in dupes:
|
|
|
399 if x[1]==y[0]:
|
|
|
400 x[0]=y[1]
|
|
|
401
|
|
|
402 for x in LH8E:
|
|
|
403 for y in dupes:
|
|
|
404 if x[1]==y[0]:
|
|
|
405 x[0]=y[1]
|
|
|
406
|
|
|
407
|
|
|
408
|
|
|
409 LH2E.sort()
|
|
|
410 LH2E=list(LH2E for LH2E,_ in itertools.groupby(LH2E))
|
|
|
411
|
|
|
412 LH8E.sort()
|
|
|
413 LH8E=list(LH8E for LH8E,_ in itertools.groupby(LH8E))
|
|
|
414
|
|
|
415 if int(per)!=-1:
|
|
|
416 percent=int(per)/100
|
|
|
417 print(percent)
|
|
|
418 print(count)
|
|
|
419
|
|
|
420 c_col_filter=round(percent*(len(LH2E[1])-2))
|
|
|
421 t_col_filter=round(percent*(len(LH8E[1])-2))
|
|
|
422
|
|
|
423 for i, _ in enumerate(LH2E):
|
|
|
424 c_cols=0
|
|
|
425 t_cols=0
|
|
|
426
|
|
|
427 c_cols=sum([1 for j in range(len(LH2E[i])-2) if int(LH2E[i][j+2])>=int(count)])
|
|
|
428 t_cols=sum([1 for j in range(len(LH8E[i])-2) if int(LH8E[i][j+2])>=int(count)])
|
|
|
429
|
|
|
430 if c_cols>=c_col_filter or t_cols>=t_col_filter:
|
|
|
431 filter_LH8E.append(LH8E[i])
|
|
|
432 filter_LH2E.append(LH2E[i])
|
|
|
433
|
|
|
434 raw_LH2E.extend(LH2E)
|
|
|
435 raw_LH8E.extend(LH8E)
|
|
|
436
|
|
|
437 ##################################################################################################################################################################################################################
|
|
|
438
|
|
|
439 def write_main(raw_LH2E, raw_LH8E, fil_LH2E, fil_LH8E, names_con, names_tre, flag, per, n1, n2):
|
|
|
440
|
|
|
441 if flag == 1 and int(per)!=-1:
|
|
|
442 fp = open('Counts/Filtered '+n2 +' Templated Counts', 'w')
|
|
|
443 fp.write("Name\t")
|
|
|
444 fp.write("Sequence")
|
|
|
445 for y in names_tre:
|
|
|
446 fp.write("\t"+y)
|
|
|
447
|
|
|
448 for x in fil_LH8E:
|
|
|
449 fp.write("\n%s" % "\t".join(x))
|
|
|
450 fp.close()
|
|
|
451
|
|
|
452 fp = open('Counts/Filtered '+n1+' Templated Counts', 'w')
|
|
|
453 fp.write("Name\t")
|
|
|
454 fp.write("Sequence")
|
|
|
455 for y in names_con:
|
|
|
456 fp.write("\t"+y)
|
|
|
457
|
|
|
458 for x in fil_LH2E:
|
|
|
459 fp.write("\n%s" % "\t".join(x))
|
|
|
460 fp.close()
|
|
|
461
|
|
|
462
|
|
|
463 if flag == 2 and int(per)!=-1:
|
|
|
464 fp = open('Counts/Filtered '+n2+' Non-Templated Counts', 'w')
|
|
|
465 fp.write("Name\t")
|
|
|
466 fp.write("Sequence")
|
|
|
467 for y in names_tre:
|
|
|
468 fp.write("\t"+y)
|
|
|
469
|
|
|
470
|
|
|
471 for x in fil_LH8E:
|
|
|
472 fp.write("\n%s" % "\t".join(x))
|
|
|
473 fp.close()
|
|
|
474
|
|
|
475 fp = open('Counts/Filtered '+n1+' Non-Templated Counts', 'w')
|
|
|
476 fp.write("Name\t")
|
|
|
477 fp.write("Sequence")
|
|
|
478 for y in names_con:
|
|
|
479 fp.write("\t"+y)
|
|
|
480
|
|
|
481 for x in fil_LH2E:
|
|
|
482 fp.write("\n%s" % "\t".join(x))
|
|
|
483 fp.close()
|
|
|
484
|
|
|
485
|
|
|
486 if flag == 1:
|
|
|
487 fp = open('Counts/Raw '+n2+' Templated Counts', 'w')
|
|
|
488 fp.write("Name\t")
|
|
|
489 fp.write("Sequence")
|
|
|
490 for y in names_tre:
|
|
|
491 fp.write("\t"+y)
|
|
|
492
|
|
|
493 for x in raw_LH8E:
|
|
|
494 fp.write("\n%s" % "\t".join(x))
|
|
|
495 fp.close()
|
|
|
496
|
|
|
497 fp = open('Counts/Raw '+n1+' Templated Counts', 'w')
|
|
|
498 fp.write("Name\t")
|
|
|
499 fp.write("Sequence")
|
|
|
500 for y in names_con:
|
|
|
501 fp.write("\t"+y)
|
|
|
502
|
|
|
503 for x in raw_LH2E:
|
|
|
504 fp.write("\n%s" % "\t".join(x))
|
|
|
505 fp.close()
|
|
|
506
|
|
|
507 if flag == 2:
|
|
|
508 fp = open('Counts/Raw '+n2+' Non-Templated Counts', 'w')
|
|
|
509 fp.write("Name\t")
|
|
|
510 fp.write("Sequence")
|
|
|
511 for y in names_tre:
|
|
|
512 fp.write("\t"+y)
|
|
|
513
|
|
|
514
|
|
|
515 for x in raw_LH8E:
|
|
|
516 fp.write("\n%s" % "\t".join(x))
|
|
|
517 fp.close()
|
|
|
518
|
|
|
519 fp = open('Counts/Raw '+n1+' Non-Templated Counts', 'w')
|
|
|
520 fp.write("Name\t")
|
|
|
521 fp.write("Sequence")
|
|
|
522 for y in names_con:
|
|
|
523 fp.write("\t"+y)
|
|
|
524
|
|
|
525 for x in raw_LH2E:
|
|
|
526 fp.write("\n%s" % "\t".join(x))
|
|
|
527 fp.close()
|
|
|
528
|
|
|
529 ####################################################################################################################################################################################################################
|
|
|
530
|
|
|
531 def ssamples(names,samp,folder,pro):
|
|
|
532
|
|
|
533 for i in range(2,len(samp[0])):
|
|
|
534
|
|
|
535 fp = open(folder+names[i-2]+'.txt','w')
|
|
|
536 fp.write("miRNA id"+"\t"+names[i-2]+"\n")
|
|
|
537
|
|
|
538 for x in samp:
|
|
|
539 fp.write("%s" % "\t".join([x[0],x[i]])+"\n")
|
|
|
540 fp.close()
|
|
|
541
|
|
|
542 ####################################################################################################################################################################################################################
|
|
|
543
|
|
|
544 def DB_write(con,name,unique_seq,sorted_uni_arms,f):
|
|
|
545
|
|
|
546 if f==1:
|
|
|
547 # Write a txt file with all the information
|
|
|
548 if con=="c":
|
|
|
549 fp = open('split1/'+name, 'w')
|
|
|
550
|
|
|
551 fp.write("%s\t%-42s\t%s\n\n" % ("Number of Reads","Name of isomir","Sequence"))
|
|
|
552 if con=="t":
|
|
|
553 fp = open('split2/'+name, 'w')
|
|
|
554 fp.write("%s\t%-42s\t%s\n\n" % ("Number of Reads","Name of isomir","Sequence"))
|
|
|
555
|
|
|
556 for i in range(len(sorted_uni_arms)):
|
|
|
557 temp = []
|
|
|
558 for j in range(len(unique_seq)):
|
|
|
559
|
|
|
560 if sorted_uni_arms[i][0] in (unique_seq[j][2].split("_")[0]+"_"+unique_seq[j][2].split("_")[1]):
|
|
|
561
|
|
|
562 temp.append(unique_seq[j])
|
|
|
563
|
|
|
564 temp = sorted(temp, key=lambda x: int(x[0].split('-')[1]), reverse=True)
|
|
|
565 fp.write("*********************************************************************************************************\n")
|
|
|
566 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]),"|"))
|
|
|
567 fp.write("*********************************************************************************************************\n\n")
|
|
|
568 [fp.write("%-8s\t%-40s\t%s\n" % (x[0].split("-")[1], x[2],x[9])) for x in temp]
|
|
|
569 fp.write("\n" + "\n")
|
|
|
570 fp.close()
|
|
|
571
|
|
|
572 if f==2:
|
|
|
573
|
|
|
574 if con=="c":
|
|
|
575 fp = open('split3/'+name, 'w')
|
|
|
576 fp.write("%s\t%-42s\t%s\n\n" % ("Number of Reads","Name of isomir","Sequence"))
|
|
|
577 if con=="t":
|
|
|
578 fp = open('split4/'+name, 'w')
|
|
|
579 fp.write("%s\t%-42s\t%s\n\n" % ("Number of Reads","Name of isomir","Sequence"))
|
|
|
580
|
|
|
581 for i in range(len(sorted_uni_arms)):
|
|
|
582 temp = []
|
|
|
583 for j in range(len(unique_seq)):
|
|
|
584 if sorted_uni_arms[i][0]==unique_seq[j][2].split("__")[0]:
|
|
|
585 temp.append(unique_seq[j])
|
|
|
586 if temp!=[]:
|
|
|
587 temp = sorted(temp, key=lambda x: int(x[0].split('-')[1]), reverse=True)
|
|
|
588 fp.write("*********************************************************************************************************\n")
|
|
|
589 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]),"|"))
|
|
|
590 fp.write("*********************************************************************************************************\n\n")
|
|
|
591 [fp.write("%-8s\t%-40s\t%s\n" % (x[0].split("-")[1], x[2],x[9])) for x in temp]
|
|
|
592 fp.write("\n" + "\n")
|
|
|
593 fp.close()
|
|
|
594
|
|
|
595
|
|
|
596 ##########################################################################################################################
|
|
|
597
|
|
|
598 def new_mat_seq(pre_unique_seq,mat_mirnas,l):
|
|
|
599
|
|
|
600 unique_iso = []
|
|
|
601 for x in pre_unique_seq:
|
|
|
602 if len(x[2].split("_"))==3:
|
|
|
603 for y in pre_unique_seq:
|
|
|
604 if x[2] in y[2] and int(x[0].split("-")[1])<int(y[0].split("-")[1]):
|
|
|
605 if any(y[2] in lst2 for lst2 in unique_iso)==False:
|
|
|
606 y[2]=">"+y[2]
|
|
|
607 unique_iso.append(y)
|
|
|
608 l.acquire()
|
|
|
609 for x in unique_iso:
|
|
|
610 mat_mirnas.append(x[2])
|
|
|
611 mat_mirnas.append(x[9])
|
|
|
612 l.release()
|
|
|
613
|
|
|
614 #########################################################################################################################
|
|
|
615
|
|
|
616 def merging_names(LH2E_copy,new):
|
|
|
617
|
|
|
618 dupes=[]
|
|
|
619 final_LH2E =[]
|
|
|
620
|
|
|
621 for num in range(len(LH2E_copy)):
|
|
|
622
|
|
|
623 if LH2E_copy[num][1] not in final_LH2E and LH2E_copy[num][0] not in final_LH2E:
|
|
|
624 final_LH2E.append(LH2E_copy[num][1])
|
|
|
625 final_LH2E.append(LH2E_copy[num][0])
|
|
|
626 else:
|
|
|
627 dupes.append(LH2E_copy[num][1])
|
|
|
628
|
|
|
629 dupes=list(set(dupes))
|
|
|
630
|
|
|
631 for i in range(len(dupes)):
|
|
|
632 dupes[i]=[dupes[i]]
|
|
|
633
|
|
|
634 for x in LH2E_copy:
|
|
|
635 for y in dupes:
|
|
|
636 if x[1]==y[0]:
|
|
|
637 fl=0
|
|
|
638 if len(y)==1:
|
|
|
639 y.append(x[0])
|
|
|
640 else:
|
|
|
641 for i in range(1,len(y)):
|
|
|
642 if y[i].split("_")[0]+"_"+y[i].split("_")[1]==x[0].split("_")[0]+"_"+x[0].split("_")[1]:
|
|
|
643 fl=1
|
|
|
644 if len(x[0])<len(y[i]):
|
|
|
645 del y[i]
|
|
|
646 y.append(x[0])
|
|
|
647 break
|
|
|
648
|
|
|
649 if fl==0:
|
|
|
650 y.append((x[0]))
|
|
|
651
|
|
|
652 for y in dupes:
|
|
|
653 if len(y)>2:
|
|
|
654 for i in range(len(y)-1,1,-1):
|
|
|
655 y[1]=y[1]+"/"+y[i]
|
|
|
656 del y[i]
|
|
|
657
|
|
|
658
|
|
|
659 for x in LH2E_copy:
|
|
|
660 for y in dupes:
|
|
|
661 if x[1]==y[0]:
|
|
|
662 x[0]=y[1]
|
|
|
663
|
|
|
664 LH2E_copy.sort()
|
|
|
665 LH2E_copy=list(LH2E_copy for LH2E_copy,_ in itertools.groupby(LH2E_copy))
|
|
|
666 new.extend(LH2E_copy)
|
|
|
667
|
|
|
668 ######################################################################################################################################################
|
|
|
669
|
|
|
670 def ssamples1(tem_names,tem_samp,non_names,non_samp,folder,pro):
|
|
|
671
|
|
|
672 for i in range(2,len(tem_samp[0])):
|
|
|
673
|
|
|
674 fp = open(folder+tem_names[i-2]+'.txt','w')
|
|
|
675 fp.write("miRNA id"+"\t"+tem_names[i-2]+"\n")
|
|
|
676
|
|
|
677 for x in tem_samp:
|
|
|
678 fp.write("%s" % "\t".join([x[0],x[i]])+"\n")
|
|
|
679
|
|
|
680 for j in range(len(non_names)):
|
|
|
681 if non_names[j]==tem_names[i-2]:
|
|
|
682 for x in non_samp:
|
|
|
683 fp.write("%s" % "\t".join([x[0],x[j+2]])+"\n")
|
|
|
684 fp.close()
|
|
|
685
|
|
|
686 #################################################################################################################################################################################################################
|
|
|
687
|
|
|
688 def download_matures(matures,org_name):
|
|
|
689
|
|
|
690 mature_mir=[]
|
|
|
691
|
|
|
692 mat_url = 'http://mirgenedb.org/fasta/'+org_name+'?mat=1'
|
|
|
693 star_url = 'http://mirgenedb.org/fasta/'+org_name+'?star=1'
|
|
|
694
|
|
|
695 data = urllib.request.urlopen(mat_url).read()
|
|
|
696 file_mirna = data.decode('utf-8')
|
|
|
697 mature_mir = file_mirna.split("\n")
|
|
|
698 mature_mir = [x.replace(">","") for x in mature_mir]
|
|
|
699 del mature_mir[-1]
|
|
|
700
|
|
|
701 data = urllib.request.urlopen(star_url).read()
|
|
|
702 file_mirna = data.decode('utf-8')
|
|
|
703 star_mir = file_mirna.split("\n")
|
|
|
704 star_mir = [x.replace(">","") for x in star_mir]
|
|
|
705 del star_mir[-1]
|
|
|
706
|
|
|
707 mature_mir.extend(star_mir)
|
|
|
708
|
|
|
709 for i in range(1,len(mature_mir),2):
|
|
|
710 mature_mir[i]=mature_mir[i].replace("U","T")
|
|
|
711
|
|
|
712 matures.extend(mature_mir)
|
|
|
713
|
|
|
714 ###################################################################################################################
|
|
|
715
|
|
|
716 def non_template_ref(sc,st,all_isoforms):
|
|
|
717
|
|
|
718 pre_uni_seq_con = list(sc)
|
|
|
719 pre_uni_seq_tre = list(st)
|
|
|
720
|
|
|
721 for x in pre_uni_seq_con:
|
|
|
722 for y in x:
|
|
|
723 if y[2] not in all_isoforms and len(y[2].split("_"))>2:
|
|
|
724 all_isoforms.append(y[2])
|
|
|
725 all_isoforms.append(y[9])
|
|
|
726
|
|
|
727 for x in pre_uni_seq_tre:
|
|
|
728 for y in x:
|
|
|
729 if y[2] not in all_isoforms and len(y[2].split("_"))>2:
|
|
|
730 all_isoforms.append(y[2])
|
|
|
731 all_isoforms.append(y[9])
|
|
|
732
|
|
|
733 ################################################################################################################################################################################################
|
|
|
734
|
|
|
735 def deseqe2(sample,mir_names,l,new_d,sample_name,sample_order):
|
|
|
736
|
|
|
737 for y in mir_names:
|
|
|
738 flag=0
|
|
|
739 for x in sample:
|
|
|
740 if y[0]==x[0]:
|
|
|
741 flag=1
|
|
|
742 break
|
|
|
743 if flag==0:
|
|
|
744 sample.append([y[0],"0",y[1]])
|
|
|
745
|
|
|
746 sample.sort(key=lambda x: x[0])
|
|
|
747 sample=list(sample for sample,_ in itertools.groupby(sample))
|
|
|
748
|
|
|
749 l.acquire()
|
|
|
750 new_d.append(sample)
|
|
|
751 sample_order.append(sample_name)
|
|
|
752 l.release()
|
|
|
753
|
|
|
754 ###############################################################################################################################################################################################
|
|
|
755
|
