changeset 12:2b09cbf138c8 draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pre_mirbase.py	Tue Oct 20 09:41:39 2020 +0000
@@ -0,0 +1,138 @@
+from itertools import groupby
+import sys
+import subprocess
+import argparse
+import time
+import urllib.request
+from multiprocessing import Process, Queue
+import itertools
+
+subprocess.call(['mkdir', 'out'])
+parser = argparse.ArgumentParser()
+
+parser.add_argument("-pos", "--positions", help="", action="store")
+parser.add_argument("-tool_dir", "--tool_directory", help="tool directory path", action="store")
+parser.add_argument("-gen", "--genome", help="tool directory path", action="store")
+parser.add_argument("-gff3", "--gff", help="",action="store")
+
+args = parser.parse_args()
+
+#=======================================================================================================================================
+
+
+#-----------------------Download and read the file hsa.gff3---------------------------------
+
+def read_url(q):
+
+
+    url = 'ftp://mirbase.org/pub/mirbase/CURRENT/genomes/'+args.gff
+    #url = 'ftp://mirbase.org/pub/mirbase/21/genomes/hsa.gff3'
+    response = urllib.request.urlopen(url)
+    data = response.read()
+    file_mirna = data.decode('utf-8')
+    file_mirna = file_mirna.split("\n")
+    q.put(file_mirna)
+
+
+def write_gff(file_mirna):
+    f = open('original_mirnas.bed', "w")
+
+    for i in range(len(file_mirna)):
+        f.write(file_mirna[i] + "\n")
+
+
+#------------------------Processed the file with mature mirnas-------------------------------
+
+
+def new_gff(file_mirna):
+
+    mirna = []   # new list with shifted mirnas
+    positions =int(args.positions)   # positions shifted
+    print(str(positions)+" positions shifted")
+    names=[]
+    # Remove lines which conatain the word "primary"
+    for i in range(len(file_mirna)):
+
+        if "primary" not in file_mirna[i]:
+            mirna.append(file_mirna[i])
+
+            if "chr" in file_mirna[i]:
+                a=file_mirna[i].split("\t")[0]
+                b=file_mirna[i].split("\t")[6]
+                c=file_mirna[i].split("=")[3].split(";")[0]
+                names.append([a,b,c])
+
+    names.sort()
+    sublists=[]
+
+    [sublists.append([item] * names.count(item)) for item in names if names.count(item)>=2]
+    sublists.sort()
+    sublists=list(sublists for sublists, _ in itertools.groupby(sublists))
+    unique_names=[[x[0][0],x[0][2]] for x in sublists]
+
+    for x in unique_names:
+        flag = 0
+        for i in range(len(mirna)):
+
+              if "chr" in mirna[i] and mirna[i].split("=")[3].split(";")[0]==x[1] and x[0]==mirna[i].split("\t")[0]:
+                 flag+=1
+                 ktr=mirna[i].split(";")[0]+";"+mirna[i].split(";")[1]+";"+mirna[i].split(";")[2]+"-{"+str(flag)+"}"+";"+mirna[i].split(";")[3]
+                 mirna[i]=ktr
+
+
+    f = open('shifted_mirnas.bed', "w")
+
+    for i in range(len(mirna)):
+
+        if "chr" in mirna[i]:
+
+            # change the name of current mirna
+            mirna_name_1 = mirna[i].split("=")[3]
+            mirna_name_2 = mirna[i].split("=")[4]
+            # mirna_name_2 = mirna_name_2.split(";")[0]
+            mirna_name_1 = mirna_name_1.split(";")[0]+"_"+mirna_name_2+"_"+mirna[i].split("\t")[0]
+            mirna[i] = mirna[i].replace("miRNA", mirna_name_1)
+
+            # Shift the position of mirna
+            start = mirna[i].split("\t")[3]
+            end = mirna[i].split("\t")[4]
+            shift_start = int(start)-positions          # shift the interval to the left
+            shift_end = int(end)+positions              # shift the interval to the right
+
+            # Replace the previous intervals with the new
+            mirna[i] = mirna[i].replace(start, str(shift_start))
+            mirna[i] = mirna[i].replace(end, str(shift_end))
+
+            f.write(mirna[i] + "\n")
+
+    f.close()
+
+#===================================================================================================================================
+
+def bedtool(genome):
+    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"])
+
+#===================================================================================================================================
+
+
+if __name__=='__main__':
+
+    starttime = time.time()
+    q = Queue()
+    p1 = Process(target=read_url(q))
+    p1.start()
+    p1.join()
+
+    file_mirna=q.get()
+
+    p2 = [Process(target=write_gff(file_mirna))]
+    p2.extend([Process(target=new_gff(file_mirna))])
+    [x.start() for x in p2]
+    [x.join() for x in p2]
+
+    p3 = Process(target=bedtool(args.genome))
+    p3.start()
+    p3.join()
+
+    print('That took {} seconds'.format(time.time() - starttime))
+