0
|
1 #!/usr/bin/env python
|
|
2
|
|
3 from galaxy import eggs
|
|
4 import sys, tempfile, os
|
|
5
|
|
6 assert sys.version_info[:2] >= (2.4)
|
|
7
|
|
8 def stop_err(msg):
|
|
9 sys.stderr.write(msg)
|
|
10 sys.exit()
|
|
11
|
|
12 def main():
|
|
13
|
|
14 out_fname = sys.argv[1]
|
|
15 in_fname = sys.argv[2]
|
|
16 chr_col = int(sys.argv[3])-1
|
|
17 coord_col = int(sys.argv[4])-1
|
|
18 track_type = sys.argv[5]
|
|
19 if track_type == 'coverage' or track_type == 'both':
|
|
20 coverage_col = int(sys.argv[6])-1
|
|
21 cname = sys.argv[7]
|
|
22 cdescription = sys.argv[8]
|
|
23 ccolor = sys.argv[9].replace('-',',')
|
|
24 cvisibility = sys.argv[10]
|
|
25 if track_type == 'snp' or track_type == 'both':
|
|
26 if track_type == 'both':
|
|
27 j = 5
|
|
28 else:
|
|
29 j = 0
|
|
30 #sname = sys.argv[7+j]
|
|
31 sdescription = sys.argv[6+j]
|
|
32 svisibility = sys.argv[7+j]
|
|
33 #ref_col = int(sys.argv[10+j])-1
|
|
34 read_col = int(sys.argv[8+j])-1
|
|
35
|
|
36
|
|
37 # Sort the input file based on chromosome (alphabetically) and start co-ordinates (numerically)
|
|
38 sorted_infile = tempfile.NamedTemporaryFile()
|
|
39 try:
|
|
40 os.system("sort -k %d,%d -k %dn -o %s %s" %(chr_col+1,chr_col+1,coord_col+1,sorted_infile.name,in_fname))
|
|
41 except Exception, exc:
|
|
42 stop_err( 'Initialization error -> %s' %str(exc) )
|
|
43
|
|
44 #generate chr list
|
|
45 sorted_infile.seek(0)
|
|
46 chr_vals = []
|
|
47 for line in file( sorted_infile.name ):
|
|
48 line = line.strip()
|
|
49 if not(line):
|
|
50 continue
|
|
51 try:
|
|
52 fields = line.split('\t')
|
|
53 chr = fields[chr_col]
|
|
54 if chr not in chr_vals:
|
|
55 chr_vals.append(chr)
|
|
56 except:
|
|
57 pass
|
|
58 if not(chr_vals):
|
|
59 stop_err("Skipped all lines as invalid.")
|
|
60
|
|
61 if track_type == 'coverage' or track_type == 'both':
|
|
62 if track_type == 'coverage':
|
|
63 fout = open( out_fname, "w" )
|
|
64 else:
|
|
65 fout = tempfile.NamedTemporaryFile()
|
|
66 fout.write('''track type=wiggle_0 name="%s" description="%s" color=%s visibility=%s\n''' \
|
|
67 % ( cname, cdescription, ccolor, cvisibility ))
|
|
68 if track_type == 'snp' or track_type == 'both':
|
|
69 fout_a = tempfile.NamedTemporaryFile()
|
|
70 fout_t = tempfile.NamedTemporaryFile()
|
|
71 fout_g = tempfile.NamedTemporaryFile()
|
|
72 fout_c = tempfile.NamedTemporaryFile()
|
|
73 fout_ref = tempfile.NamedTemporaryFile()
|
|
74
|
|
75 fout_a.write('''track type=wiggle_0 name="%s" description="%s" color=%s visibility=%s\n''' \
|
|
76 % ( "Track A", sdescription, '255,0,0', svisibility ))
|
|
77 fout_t.write('''track type=wiggle_0 name="%s" description="%s" color=%s visibility=%s\n''' \
|
|
78 % ( "Track T", sdescription, '0,255,0', svisibility ))
|
|
79 fout_g.write('''track type=wiggle_0 name="%s" description="%s" color=%s visibility=%s\n''' \
|
|
80 % ( "Track G", sdescription, '0,0,255', svisibility ))
|
|
81 fout_c.write('''track type=wiggle_0 name="%s" description="%s" color=%s visibility=%s\n''' \
|
|
82 % ( "Track C", sdescription, '255,0,255', svisibility ))
|
|
83
|
|
84
|
|
85 sorted_infile.seek(0)
|
|
86 for line in file( sorted_infile.name ):
|
|
87 line = line.strip()
|
|
88 if not(line):
|
|
89 continue
|
|
90 try:
|
|
91 fields = line.split('\t')
|
|
92 chr = fields[chr_col]
|
|
93 start = int(fields[coord_col])
|
|
94 assert start > 0
|
|
95 except:
|
|
96 continue
|
|
97 try:
|
|
98 ind = chr_vals.index(chr) #encountered chr for the 1st time
|
|
99 del chr_vals[ind]
|
|
100 prev_start = ''
|
|
101 header = "variableStep chrom=%s\n" %(chr)
|
|
102 if track_type == 'coverage' or track_type == 'both':
|
|
103 coverage = int(fields[coverage_col])
|
|
104 line1 = "%s\t%s\n" %(start,coverage)
|
|
105 fout.write("%s%s" %(header,line1))
|
|
106 if track_type == 'snp' or track_type == 'both':
|
|
107 a = t = g = c = 0
|
|
108 fout_a.write("%s" %(header))
|
|
109 fout_t.write("%s" %(header))
|
|
110 fout_g.write("%s" %(header))
|
|
111 fout_c.write("%s" %(header))
|
|
112 try:
|
|
113 #ref_nt = fields[ref_col].capitalize()
|
|
114 read_nt = fields[read_col].capitalize()
|
|
115 try:
|
|
116 nt_ind = ['A','T','G','C'].index(read_nt)
|
|
117 if nt_ind == 0:
|
|
118 a+=1
|
|
119 elif nt_ind == 1:
|
|
120 t+=1
|
|
121 elif nt_ind == 2:
|
|
122 g+=1
|
|
123 else:
|
|
124 c+=1
|
|
125 except ValueError:
|
|
126 pass
|
|
127 except:
|
|
128 pass
|
|
129 prev_start = start
|
|
130 except ValueError:
|
|
131 if start != prev_start:
|
|
132 if track_type == 'coverage' or track_type == 'both':
|
|
133 coverage = int(fields[coverage_col])
|
|
134 fout.write("%s\t%s\n" %(start,coverage))
|
|
135 if track_type == 'snp' or track_type == 'both':
|
|
136 if a:
|
|
137 fout_a.write("%s\t%s\n" %(prev_start,a))
|
|
138 if t:
|
|
139 fout_t.write("%s\t%s\n" %(prev_start,t))
|
|
140 if g:
|
|
141 fout_g.write("%s\t%s\n" %(prev_start,g))
|
|
142 if c:
|
|
143 fout_c.write("%s\t%s\n" %(prev_start,c))
|
|
144 a = t = g = c = 0
|
|
145 try:
|
|
146 #ref_nt = fields[ref_col].capitalize()
|
|
147 read_nt = fields[read_col].capitalize()
|
|
148 try:
|
|
149 nt_ind = ['A','T','G','C'].index(read_nt)
|
|
150 if nt_ind == 0:
|
|
151 a+=1
|
|
152 elif nt_ind == 1:
|
|
153 t+=1
|
|
154 elif nt_ind == 2:
|
|
155 g+=1
|
|
156 else:
|
|
157 c+=1
|
|
158 except ValueError:
|
|
159 pass
|
|
160 except:
|
|
161 pass
|
|
162 prev_start = start
|
|
163 else:
|
|
164 if track_type == 'snp' or track_type == 'both':
|
|
165 try:
|
|
166 #ref_nt = fields[ref_col].capitalize()
|
|
167 read_nt = fields[read_col].capitalize()
|
|
168 try:
|
|
169 nt_ind = ['A','T','G','C'].index(read_nt)
|
|
170 if nt_ind == 0:
|
|
171 a+=1
|
|
172 elif nt_ind == 1:
|
|
173 t+=1
|
|
174 elif nt_ind == 2:
|
|
175 g+=1
|
|
176 else:
|
|
177 c+=1
|
|
178 except ValueError:
|
|
179 pass
|
|
180 except:
|
|
181 pass
|
|
182
|
|
183 if track_type == 'snp' or track_type == 'both':
|
|
184 if a:
|
|
185 fout_a.write("%s\t%s\n" %(prev_start,a))
|
|
186 if t:
|
|
187 fout_t.write("%s\t%s\n" %(prev_start,t))
|
|
188 if g:
|
|
189 fout_g.write("%s\t%s\n" %(prev_start,g))
|
|
190 if c:
|
|
191 fout_c.write("%s\t%s\n" %(prev_start,c))
|
|
192
|
|
193 fout_a.seek(0)
|
|
194 fout_g.seek(0)
|
|
195 fout_t.seek(0)
|
|
196 fout_c.seek(0)
|
|
197
|
|
198 if track_type == 'snp':
|
|
199 os.system("cat %s %s %s %s >> %s" %(fout_a.name,fout_t.name,fout_g.name,fout_c.name,out_fname))
|
|
200 elif track_type == 'both':
|
|
201 fout.seek(0)
|
|
202 os.system("cat %s %s %s %s %s | cat > %s" %(fout.name,fout_a.name,fout_t.name,fout_g.name,fout_c.name,out_fname))
|
|
203 if __name__ == "__main__":
|
|
204 main() |