|
0
|
1 """
|
|
|
2 Converts any type of FASTQ file to Sanger type and makes small adjustments if necessary.
|
|
|
3
|
|
|
4 usage: %prog [options]
|
|
|
5 -i, --input=i: Input FASTQ candidate file
|
|
|
6 -r, --origType=r: Original type
|
|
|
7 -a, --allOrNot=a: Whether or not to check all blocks
|
|
|
8 -b, --blocks=b: Number of blocks to check
|
|
|
9 -o, --output=o: Output file
|
|
|
10
|
|
|
11 usage: %prog input_file oroutput_file
|
|
|
12 """
|
|
|
13
|
|
|
14 import math, sys
|
|
|
15 from galaxy import eggs
|
|
|
16 import pkg_resources; pkg_resources.require( "bx-python" )
|
|
|
17 from bx.cookbook import doc_optparse
|
|
|
18
|
|
|
19 def stop_err( msg ):
|
|
|
20 sys.stderr.write( "%s\n" % msg )
|
|
|
21 sys.exit()
|
|
|
22
|
|
|
23 def all_bases_valid(seq):
|
|
|
24 """Confirm that the sequence contains only bases"""
|
|
|
25 valid_bases = ['a', 'A', 'c', 'C', 'g', 'G', 't', 'T', 'N']
|
|
|
26 for base in seq:
|
|
|
27 if base not in valid_bases:
|
|
|
28 return False
|
|
|
29 return True
|
|
|
30
|
|
|
31 def __main__():
|
|
|
32 #Parse Command Line
|
|
|
33 options, args = doc_optparse.parse( __doc__ )
|
|
|
34 orig_type = options.origType
|
|
|
35 if orig_type == 'sanger' and options.allOrNot == 'not':
|
|
|
36 max_blocks = int(options.blocks)
|
|
|
37 else:
|
|
|
38 max_blocks = -1
|
|
|
39 fin = file(options.input, 'r')
|
|
|
40 fout = file(options.output, 'w')
|
|
|
41 range_min = 1000
|
|
|
42 range_max = -5
|
|
|
43 block_num = 0
|
|
|
44 bad_blocks = 0
|
|
|
45 base_len = -1
|
|
|
46 line_count = 0
|
|
|
47 lines = []
|
|
|
48 line = fin.readline()
|
|
|
49 while line:
|
|
|
50 if line.strip() and max_blocks >= 0 and block_num > 0 and orig_type == 'sanger' and block_num >= max_blocks:
|
|
|
51 fout.write(line)
|
|
|
52 if line_count % 4 == 0:
|
|
|
53 block_num += 1
|
|
|
54 line_count += 1
|
|
|
55 elif line.strip():
|
|
|
56 # the line that starts a block, with a name
|
|
|
57 if line_count % 4 == 0 and line.startswith('@'):
|
|
|
58 lines.append(line)
|
|
|
59 else:
|
|
|
60 # if we expect a sequence of bases
|
|
|
61 if line_count % 4 == 1 and all_bases_valid(line.strip()):
|
|
|
62 lines.append(line)
|
|
|
63 base_len = len(line.strip())
|
|
|
64 # if we expect the second name line
|
|
|
65 elif line_count % 4 == 2 and line.startswith('+'):
|
|
|
66 lines.append(line)
|
|
|
67 # if we expect a sequence of qualities and it's the expected length
|
|
|
68 elif line_count % 4 == 3:
|
|
|
69 split_line = line.strip().split()
|
|
|
70 # decimal qualities
|
|
|
71 if len(split_line) == base_len:
|
|
|
72 # convert
|
|
|
73 phred_list = []
|
|
|
74 for ch in split_line:
|
|
|
75 int_ch = int(ch)
|
|
|
76 if int_ch < range_min:
|
|
|
77 range_min = int_ch
|
|
|
78 if int_ch > range_max:
|
|
|
79 range_max = int_ch
|
|
|
80 if int_ch >= 0 and int_ch <= 93:
|
|
|
81 phred_list.append(chr(int_ch + 33))
|
|
|
82 # make sure we haven't lost any quality values
|
|
|
83 if len(phred_list) == base_len:
|
|
|
84 # print first three lines
|
|
|
85 for l in lines:
|
|
|
86 fout.write(l)
|
|
|
87 # print converted quality line
|
|
|
88 fout.write(''.join(phred_list))
|
|
|
89 # reset
|
|
|
90 lines = []
|
|
|
91 base_len = -1
|
|
|
92 # abort if so
|
|
|
93 else:
|
|
|
94 bad_blocks += 1
|
|
|
95 lines = []
|
|
|
96 base_len = -1
|
|
|
97 # ascii qualities
|
|
|
98 elif len(split_line[0]) == base_len:
|
|
|
99 qualities = []
|
|
|
100 # print converted quality line
|
|
|
101 if orig_type == 'illumina':
|
|
|
102 for c in line.strip():
|
|
|
103 if ord(c) - 64 < range_min:
|
|
|
104 range_min = ord(c) - 64
|
|
|
105 if ord(c) - 64 > range_max:
|
|
|
106 range_max = ord(c) - 64
|
|
|
107 if ord(c) < 64 or ord(c) > 126:
|
|
|
108 bad_blocks += 1
|
|
|
109 base_len = -1
|
|
|
110 lines = []
|
|
|
111 break
|
|
|
112 else:
|
|
|
113 qualities.append( chr( ord(c) - 31 ) )
|
|
|
114 quals = ''.join(qualities)
|
|
|
115 elif orig_type == 'solexa':
|
|
|
116 for c in line.strip():
|
|
|
117 if ord(c) - 64 < range_min:
|
|
|
118 range_min = ord(c) - 64
|
|
|
119 if ord(c) - 64 > range_max:
|
|
|
120 range_max = ord(c) - 64
|
|
|
121 if ord(c) < 59 or ord(c) > 126:
|
|
|
122 bad_blocks += 1
|
|
|
123 base_len = -1
|
|
|
124 lines = []
|
|
|
125 break
|
|
|
126 else:
|
|
|
127 p = 10.0**( ( ord(c) - 64 ) / -10.0 ) / ( 1 + 10.0**( ( ord(c) - 64 ) / -10.0 ) )
|
|
|
128 qualities.append( chr( int( -10.0*math.log10( p ) ) + 33 ) )
|
|
|
129 quals = ''.join(qualities)
|
|
|
130 else: # 'sanger'
|
|
|
131 for c in line.strip():
|
|
|
132 if ord(c) - 33 < range_min:
|
|
|
133 range_min = ord(c) - 33
|
|
|
134 if ord(c) - 33 > range_max:
|
|
|
135 range_max = ord(c) - 33
|
|
|
136 if ord(c) < 33 or ord(c) > 126:
|
|
|
137 bad_blocks += 1
|
|
|
138 base_len = -1
|
|
|
139 lines = []
|
|
|
140 break
|
|
|
141 else:
|
|
|
142 qualities.append(c)
|
|
|
143 quals = ''.join(qualities)
|
|
|
144 # make sure we don't have bad qualities
|
|
|
145 if len(quals) == base_len:
|
|
|
146 # print first three lines
|
|
|
147 for l in lines:
|
|
|
148 fout.write(l)
|
|
|
149 # print out quality line
|
|
|
150 fout.write(quals+'\n')
|
|
|
151 # reset
|
|
|
152 lines = []
|
|
|
153 base_len = -1
|
|
|
154 else:
|
|
|
155 bad_blocks += 1
|
|
|
156 base_len = -1
|
|
|
157 lines = []
|
|
|
158 # mark the successful end of a block
|
|
|
159 block_num += 1
|
|
|
160 line_count += 1
|
|
|
161 line = fin.readline()
|
|
|
162 fout.close()
|
|
|
163 fin.close()
|
|
|
164 if range_min != 1000 and range_min != -5:
|
|
|
165 outmsg = 'The range of quality values found were: %s to %s' % (range_min, range_max)
|
|
|
166 else:
|
|
|
167 outmsg = ''
|
|
|
168 if bad_blocks > 0:
|
|
|
169 outmsg += '\nThere were %s bad blocks skipped' % (bad_blocks)
|
|
|
170 sys.stdout.write(outmsg)
|
|
|
171
|
|
|
172 if __name__=="__main__": __main__() |
