/tools/next_gen_conversion/fastq_gen_conv.py

  1"""
  2Converts any type of FASTQ file to Sanger type  and makes small adjustments if necessary.
  3
  4usage: %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
 11usage: %prog input_file oroutput_file
 12"""
 13
 14import math, sys
 15from galaxy import eggs
 16import pkg_resources; pkg_resources.require( "bx-python" )
 17from bx.cookbook import doc_optparse
 18
 19def stop_err( msg ):
 20    sys.stderr.write( "%s\n" % msg )
 21    sys.exit()
 22    
 23def 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
 31def __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
172if __name__=="__main__": __main__()