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