annotate armdb_mirbase.py @ 4:fbb1f770ab77 draft

Uploaded
author glogobyte
date Wed, 13 Oct 2021 16:23:59 +0000
parents 6c267b9256fa
children 40c50594116c
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
1 import subprocess
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
2 import argparse
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
3 import time
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
4 import urllib.request
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
5 from multiprocessing import Process, Queue
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
6 import itertools
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
7
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
8
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
9
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
10 #---------------------------------Arguments------------------------------------------
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
11
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
12 subprocess.call(['mkdir', 'out'])
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
13 parser = argparse.ArgumentParser()
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
14
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
15 parser.add_argument("-pos", "--positions", help="", action="store")
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
16 parser.add_argument("-tool_dir", "--tool_directory", help="tool directory path", action="store")
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
17 parser.add_argument("-gen", "--genome", help="tool directory path", action="store")
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
18 parser.add_argument("-gff3", "--gff", help="",action="store")
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
19 args = parser.parse_args()
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
20
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
21 #-----------------------Download and read of the gff3 file---------------------------------
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
22
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
23 def read_url(q):
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
24
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
25 url = 'ftp://mirbase.org/pub/mirbase/CURRENT/genomes/'+args.gff
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
26 data = urllib.request.urlopen(url).read()
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
27 file_mirna = data.decode('utf-8')
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
28 file_mirna = file_mirna.split("\n")
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
29 q.put(file_mirna)
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
30
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
31
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
32 #-----------------------Export of the original gff3 file---------------------------------
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
33
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
34 def write_gff(file_mirna):
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
35 f = open('original_mirnas.bed', "w")
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
36
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
37 for i in range(len(file_mirna)):
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
38 f.write(file_mirna[i] + "\n")
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
39
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
40 #------------------------Process and export of the file with mature mirnas-------------------------------
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
41
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
42 def new_gff(file_mirna):
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
43
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
44 mirna = [] # new list with shifted mirnas
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
45 positions =int(args.positions) # positions shifted
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
46 print(str(positions)+" positions shifted")
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
47 names=[]
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
48
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
49 for i in range(len(file_mirna)):
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
50
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
51 # Remove lines which conatain the word "primary"
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
52 if "primary" not in file_mirna[i]:
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
53 mirna.append(file_mirna[i])
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
54
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
55 # Check if the line starts with "chr"
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
56 if "chr" in file_mirna[i]:
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
57 a=file_mirna[i].split("\t")[0]
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
58 b=file_mirna[i].split("\t")[6]
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
59 c=file_mirna[i].split("=")[3].split(";")[0]
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
60 names.append([a,b,c])
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
61
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
62 names.sort()
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
63 sublists=[]
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
64
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
65 [sublists.append([item] * names.count(item)) for item in names if names.count(item)>=2]
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
66 sublists.sort()
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
67 sublists=list(sublists for sublists, _ in itertools.groupby(sublists))
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
68 unique_names=[[x[0][0],x[0][2]] for x in sublists]
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
69
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
70 for x in unique_names:
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
71 flag = 0
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
72 for i in range(len(mirna)):
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
73
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
74 if "chr" in mirna[i] and mirna[i].split("=")[3].split(";")[0]==x[1] and x[0]==mirna[i].split("\t")[0]:
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
75 flag+=1
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
76 ktr=mirna[i].split(";")[0]+";"+mirna[i].split(";")[1]+";"+mirna[i].split(";")[2]+"-{"+str(flag)+"}"+";"+mirna[i].split(";")[3]
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
77 mirna[i]=ktr
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
78
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
79
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
80 f = open('shifted_mirnas.bed', "w")
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
81
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
82 for i in range(len(mirna)):
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
83
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
84 if "chr" in mirna[i]:
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
85
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
86 # change the name of current mirna
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
87 mirna_name_1 = mirna[i].split("=")[3]
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
88 mirna_name_2 = mirna[i].split("=")[4]
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
89 mirna_name_1 = mirna_name_1.split(";")[0]+"_"+mirna_name_2+"_"+mirna[i].split("\t")[0]
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
90 mirna[i] = mirna[i].replace("miRNA", mirna_name_1)
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
91
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
92 # Shift the position of mirna
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
93 start = mirna[i].split("\t")[3]
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
94 end = mirna[i].split("\t")[4]
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
95 shift_start = int(start)-positions # shift the interval to the left
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
96 shift_end = int(end)+positions # shift the interval to the right
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
97
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
98 # Replace the previous intervals with the new
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
99 mirna[i] = mirna[i].replace(start, str(shift_start))
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
100 mirna[i] = mirna[i].replace(end, str(shift_end))
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
101
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
102 f.write(mirna[i] + "\n")
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
103
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
104 f.close()
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
105
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
106
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
107 #------------------------Extract the sequences of the Custom Arms with getfasta tool-------------------------------
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
108
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
109 def bedtool(genome):
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
110 subprocess.call(["bedtools", "getfasta", "-fi", "/cvmfs/data.galaxyproject.org/byhand/"+genome+"/sam_index/"+genome+".fa", "-bed", "shifted_mirnas.bed", "-fo", "new_ref.fa", "-name", "-s"])
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
111 #===================================================================================================================================
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
112
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
113
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
114 if __name__=='__main__':
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
115
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
116 starttime = time.time()
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
117 q = Queue()
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
118
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
119 # Read the original gff3 file
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
120 p1 = Process(target=read_url(q))
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
121 p1.start()
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
122 p1.join()
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
123 file_mirna=q.get()
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
124
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
125 # Export the original gff3 file
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
126 p2 = [Process(target=write_gff(file_mirna))]
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
127
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
128 # Create the new gff3 file
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
129 p2.extend([Process(target=new_gff(file_mirna))])
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
130 [x.start() for x in p2]
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
131 [x.join() for x in p2]
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
132
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
133 # Extract the sequences of the Custom Arms
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
134 p3 = Process(target=bedtool(args.genome))
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
135 p3.start()
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
136 p3.join()
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
137
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
138 print('The total runtime was {} seconds'.format(time.time() - starttime))
6c267b9256fa Uploaded
glogobyte
parents:
diff changeset
139