Mercurial > repos > devteam > fastqsolexa_to_fasta_qual

comparison fastqsolexa_to_fasta_qual.py @ 0:b471da1d7303 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Imported from capsule None
author devteam
date Mon, 19 May 2014 10:59:35 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:b471da1d7303
1 #!/usr/bin/env python
2
3 """
4 convert fastqsolexa file to separated sequence and quality files.
5
6 assume each sequence and quality score are contained in one line
7 the order should be:
8 1st line: @title_of_seq
9 2nd line: nucleotides
10 3rd line: +title_of_qualityscore (might be skipped)
11 4th line: quality scores
12 (in three forms: a. digits, b. ASCII codes, the first char as the coding base, c. ASCII codes without the first char.)
13
14 Usage:
15 %python fastqsolexa_to_fasta_qual.py <your_fastqsolexa_filename> <output_seq_filename> <output_score_filename>
16 """
17
18 import sys, os
19 from math import *
20
21 assert sys.version_info[:2] >= ( 2, 4 )
22
23 def stop_err( msg ):
24 sys.stderr.write( "%s" % msg )
25 sys.exit()
26
27 def __main__():
28 infile_name = sys.argv[1]
29 outfile_seq = open( sys.argv[2], 'w' )
30 outfile_score = open( sys.argv[3], 'w' )
31 datatype = sys.argv[4]
32 seq_title_startswith = ''
33 qual_title_startswith = ''
34 default_coding_value = 64
35 fastq_block_lines = 0
36
37 for i, line in enumerate( file( infile_name ) ):
38 line = line.rstrip()
39 if not line or line.startswith( '#' ):
40 continue
41 fastq_block_lines = ( fastq_block_lines + 1 ) % 4
42 line_startswith = line[0:1]
43 if fastq_block_lines == 1:
44 # first line is @title_of_seq
45 if not seq_title_startswith:
46 seq_title_startswith = line_startswith
47 if line_startswith != seq_title_startswith:
48 outfile_seq.close()
49 outfile_score.close()
50 stop_err( 'Invalid fastqsolexa format at line %d: %s.' % ( i + 1, line ) )
51 read_title = line[1:]
52 outfile_seq.write( '>%s\n' % line[1:] )
53 elif fastq_block_lines == 2:
54 # second line is nucleotides
55 read_length = len( line )
56 outfile_seq.write( '%s\n' % line )
57 elif fastq_block_lines == 3:
58 # third line is +title_of_qualityscore ( might be skipped )
59 if not qual_title_startswith:
60 qual_title_startswith = line_startswith
61 if line_startswith != qual_title_startswith:
62 outfile_seq.close()
63 outfile_score.close()
64 stop_err( 'Invalid fastqsolexa format at line %d: %s.' % ( i + 1, line ) )
65 quality_title = line[1:]
66 if quality_title and read_title != quality_title:
67 outfile_seq.close()
68 outfile_score.close()
69 stop_err( 'Invalid fastqsolexa format at line %d: sequence title "%s" differes from score title "%s".' % ( i + 1, read_title, quality_title ) )
70 if not quality_title:
71 outfile_score.write( '>%s\n' % read_title )
72 else:
73 outfile_score.write( '>%s\n' % line[1:] )
74 else:
75 # fourth line is quality scores
76 qual = ''
77 fastq_integer = True
78 # peek: ascii or digits?
79 val = line.split()[0]
80 try:
81 check = int( val )
82 fastq_integer = True
83 except:
84 fastq_integer = False
85
86 if fastq_integer:
87 # digits
88 qual = line
89 else:
90 # ascii
91 quality_score_length = len( line )
92 if quality_score_length == read_length + 1:
93 # first char is qual_score_startswith
94 qual_score_startswith = ord( line[0:1] )
95 line = line[1:]
96 elif quality_score_length == read_length:
97 qual_score_startswith = default_coding_value
98 else:
99 stop_err( 'Invalid fastqsolexa format at line %d: the number of quality scores ( %d ) is not the same as bases ( %d ).' % ( i + 1, quality_score_length, read_length ) )
100 for j, char in enumerate( line ):
101 score = ord( char ) - qual_score_startswith # 64
102 qual = "%s%s " % ( qual, str( score ) )
103 outfile_score.write( '%s\n' % qual )
104
105 outfile_seq.close()
106 outfile_score.close()
107
108 if __name__ == "__main__": __main__()
109