|
56
|
1 #!/usr/bin/env python
|
|
|
2 # -*- coding: utf-8 -*-
|
|
|
3
|
|
|
4 import sys
|
|
|
5 import os
|
|
|
6 from read_file import *
|
|
|
7 from read_s_file import *
|
|
|
8 import random
|
|
|
9 import string
|
|
|
10
|
|
|
11 type_input = sys.argv[1]
|
|
|
12 seq_file = sys.argv[2]
|
|
|
13 ref_file = sys.argv[3]
|
|
|
14 shift = sys.argv[4]
|
|
|
15 length = sys.argv[5]
|
|
|
16 t_end = sys.argv[6]
|
|
|
17 map_type = sys.argv[7]
|
|
|
18 output_file = sys.argv[8]
|
|
|
19
|
|
|
20
|
|
|
21 if map_type!="default":
|
|
|
22 s = ""
|
|
|
23 s = s+"-v "+sys.argv[9]
|
|
|
24 s = s+" -5 "+sys.argv[10]
|
|
|
25 s = s+" -3 "+sys.argv[11]
|
|
|
26 s = s+" -k "+sys.argv[12]
|
|
|
27 if sys.argv[13]:
|
|
|
28 s = s+" -a"
|
|
|
29 if int(sys.argv[14])>=1:
|
|
|
30 s = s+" -m "+sys.argv[14]
|
|
|
31 if sys.argv[15]:
|
|
|
32 s = s+" --best --strata "
|
|
|
33
|
|
|
34 else:
|
|
|
35 s = "-v 3 -a --best --strata "
|
|
|
36
|
|
|
37 ospath = os.path.realpath(sys.argv[0])
|
|
|
38 ost = ospath.split('/')
|
|
|
39 syspath = ""
|
|
|
40 for i in range(len(ost)-1):
|
|
|
41 syspath = syspath+ost[i].strip()
|
|
|
42 syspath = syspath+'/'
|
|
|
43
|
|
|
44 rs = ''.join(random.sample(string.ascii_letters + string.digits, 8))
|
|
|
45 os.system("mkdir "+syspath+"output_"+rs)
|
|
|
46
|
|
|
47 syspathrs = syspath+"output_"+rs+'/'
|
|
|
48
|
|
|
49 os.system("bowtie-build -f "+ref_file+" "+syspathrs+"ref > "+syspathrs+"log.txt")
|
|
|
50
|
|
|
51 os.system("cp "+seq_file+" "+syspathrs+"seq0.fa")
|
|
|
52
|
|
|
53 if type_input == "fasta":
|
|
|
54 tp = 'fasta'
|
|
|
55 if type_input == "fastq":
|
|
|
56 tp = 'fastq'
|
|
|
57
|
|
|
58 k = 0
|
|
|
59 print(type_input)
|
|
|
60 while(True):
|
|
|
61 if type_input == "fasta":
|
|
|
62 os.system("bowtie "+s+"-f "+syspathrs+"ref"+" "+syspathrs+"seq"+str(k)+".fa --quiet -S > "+syspathrs+"map"+str(k)+".sam")
|
|
|
63 if type_input == "fastq":
|
|
|
64 os.system("bowtie "+s+"-q "+syspathrs+"ref"+" "+syspathrs+"seq"+str(k)+".fa --quiet -S > "+syspathrs+"map"+str(k)+".sam")
|
|
|
65 os.system("samtools view -Sb -F 0xfff "+syspathrs+"map"+str(k)+".sam > "+syspathrs+"mapped"+str(k)+".bam 2>"+syspathrs+"log.txt") #get mapped reads
|
|
|
66 os.system("samtools view -Sb -f 0x4 "+syspathrs+"map"+str(k)+".sam > "+syspathrs+"umapped"+str(k)+".bam 2>"+syspathrs+"log.txt") #get unmapped reads
|
|
|
67 os.system("samtools view -Sb -f 0x10 "+syspathrs+"map"+str(k)+".sam > "+syspathrs+"rmapped"+str(k)+".bam 2>"+syspathrs+"log.txt") #get reversed mapped reads
|
|
|
68 os.system("samtools merge -f "+syspathrs+"unmapped"+str(k)+".bam "+syspathrs+"umapped"+str(k)+".bam "+syspathrs+"rmapped"+str(k)+".bam") #get reversed mapped reads
|
|
|
69 os.system("samtools view -h -o "+syspathrs+"unmapped"+str(k)+".sam "+syspathrs+"unmapped"+str(k)+".bam") #get reversed mapped reads
|
|
|
70 if k>0:
|
|
|
71 os.system("samtools view -h -o "+syspathrs+"mapped"+str(k)+".sam "+syspathrs+"mapped"+str(k)+".bam") #get reversed mapped reads
|
|
|
72 os.system("cut -f 1 "+syspathrs+"unmapped"+str(k)+".sam > "+syspathrs+"unmapped"+str(k)+".txt")
|
|
|
73 os.system("cut -f 1 "+syspathrs+"mapped"+str(k)+".sam > "+syspathrs+"mapped"+str(k)+".txt")
|
|
|
74 os.system("python "+syspath+"remove_map.py "+syspathrs+"unmapped"+str(k)+".txt "+syspathrs+"mapped"+str(k)+".txt "+syspathrs+"runmapped"+str(k)+".txt")
|
|
|
75 os.system("rm "+syspathrs+"mapped"+str(k)+".sam")
|
|
|
76 os.system("rm "+syspathrs+"mapped"+str(k)+".txt")
|
|
|
77 os.system("rm "+syspathrs+"unmapped"+str(k)+".txt")
|
|
|
78 else:
|
|
|
79 os.system("cut -f 1 "+syspathrs+"unmapped"+str(k)+".sam > "+syspathrs+"runmapped"+str(k)+".txt")
|
|
|
80
|
|
|
81 os.system("rm "+syspathrs+"unmapped"+str(k)+".bam")
|
|
|
82 os.system("rm "+syspathrs+"umapped"+str(k)+".bam")
|
|
|
83 os.system("rm "+syspathrs+"rmapped"+str(k)+".bam")
|
|
|
84 os.system("python "+syspath+"seq_track.py "+syspathrs+"runmapped"+str(k)+".txt "+syspathrs+"seq"+str(k)+".fa "+syspathrs+"unmap_seq"+str(k)+".fa "+tp) #get unmapped sequence
|
|
|
85 os.system("python "+syspath+"truncate.py "+syspathrs+"unmap_seq"+str(k)+".fa "+shift+" "+syspathrs+"seq"+str(k+1)+".fa "+length+" "+t_end) #truncate unmapped sequence
|
|
|
86 os.system("rm "+syspathrs+"seq"+str(k)+".fa") #Remove sequences being mapped
|
|
|
87 os.system("rm "+syspathrs+"map"+str(k)+".sam") #Remove mapping file
|
|
|
88 os.system("rm "+syspathrs+"unmap_seq"+str(k)+".fa") #Remove unmapped sequnce
|
|
|
89 os.system("rm "+syspathrs+"runmapped"+str(k)+".txt")
|
|
|
90 os.system("rm "+syspathrs+"unmapped"+str(k)+".sam")
|
|
|
91
|
|
|
92 os.system("wc -l "+syspathrs+"seq"+str(k+1)+".fa > "+syspathrs+"count"+str(k+1)+".txt")
|
|
|
93 c = read_sp_file(syspathrs+"count"+str(k+1)+".txt")
|
|
|
94 if c[0][0] == '0': #If no reads is in the sequence file, stop
|
|
|
95 os.system("rm "+syspathrs+"count"+str(k+1)+".txt")
|
|
|
96 os.system("rm "+syspathrs+"seq"+str(k+1)+".fa")
|
|
|
97 break
|
|
|
98 os.system("rm "+syspathrs+"count"+str(k+1)+".txt")
|
|
|
99 k = k+1
|
|
|
100
|
|
|
101 ss = ""
|
|
|
102 for i in range(0,k+1):
|
|
|
103 ss = ss+" "+syspathrs+"mapped"+str(i)+".bam"
|
|
|
104
|
|
|
105
|
|
|
106 os.system("samtools merge -f "+output_file+" "+ss)
|
|
|
107 #print("samtools merge mapped_all.bam"+ss)
|
|
|
108 os.system("rm "+syspathrs+"mapped*.bam")
|
|
|
109 os.system("rm "+syspathrs+"ref*")
|
|
|
110 os.system("rm -r "+syspathrs)
|
|
|
111
|
|
|
112
|
