#NEWLINE# def main( argv = None ) : #NEWLINE# #INDENT# 'script main.\n\n    parses command line options in sys.argv, unless *argv* is given.\n    ' #NEWLINE# if ( not argv ) : #NEWLINE# #INDENT# argv = sys . argv #NEWLINE# #UNINDENT# parser = E . OptionParser ( version = '%prog version: $Id$' , usage = globals ( ) [ '__doc__' ] ) #NEWLINE# parser . add_option ( '-m' , '--method' , dest = 'method' , type = 'choice' , choices = ( 'apply' , 'change-format' , 'renumber-reads' , 'sample' , 'sort' , 'trim3' , 'trim5' , 'unique' , 'reverse-complement' , 'grep' ) , help = 'method to apply [%default]' ) #NEWLINE# parser . add_option ( '--target-format' , dest = 'target_format' , type = 'choice' , choices = ( 'sanger' , 'solexa' , 'phred64' , 'integer' , 'illumina-1.8' ) , help = 'guess quality score format and set quality scores to format [default=%default].' ) #NEWLINE# parser . add_option ( '--guess-format' , dest = 'guess_format' , type = 'choice' , choices = ( 'sanger' , 'solexa' , 'phred64' , 'integer' , 'illumina-1.8' ) , help = 'quality score format to assume if ambiguous [default=%default].' ) #NEWLINE# parser . add_option ( '--sample-size' , dest = 'sample_size' , type = 'float' , help = 'proportion of reads to sample. Provide a proportion of reads to sample, e.g. 0.1 for 10%, 0.5 for 50%, etc [default=%default].' ) #NEWLINE# parser . add_option ( '--pair-fastq-file' , dest = 'pair' , type = 'string' , help = 'if data is paired, filename with second pair. Implemented for sampling [default=%default].' ) #NEWLINE# parser . add_option ( '--map-tsv-file' , dest = 'map_tsv_file' , type = 'string' , help = 'filename with tab-separated identifiers mapping for method apply [default=%default].' ) #NEWLINE# parser . add_option ( '--num-bases' , dest = 'nbases' , type = 'int' , help = 'number of bases to trim [default=%default].' ) #NEWLINE# parser . add_option ( '--seed' , dest = 'seed' , type = 'int' , help = 'seed for random number generator [default=%default].' ) #NEWLINE# parser . add_option ( '--pattern-identifier' , dest = 'renumber_pattern' , type = 'string' , help = 'rename reads in file by pattern [default=%default]' ) #NEWLINE# parser . add_option ( '--grep-pattern' , dest = 'grep_pattern' , type = 'string' , help = 'subset to reads matching pattern [default=%default]' ) #NEWLINE# parser . set_defaults ( method = None , change_format = None , guess_format = None , sample_size = 0.1 , nbases = 0 , pair = None , apply = None , seed = None , renumber_pattern = 'read_%010i' , grep_pattern = '.*' ) #NEWLINE# ( options , args ) = E . Start ( parser , argv = argv , add_output_options = True ) #NEWLINE# c = E . Counter ( ) #NEWLINE# if ( options . method == 'change-format' ) : #NEWLINE# #INDENT# for record in Fastq . iterate_convert ( options . stdin , format = options . target_format , guess = options . guess_format ) : #NEWLINE# #INDENT# c . input += 1 #NEWLINE# options . stdout . write ( ( '%s\n' % record ) ) #NEWLINE# c . output += 1 #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'grep' ) : #NEWLINE# #INDENT# for record in Fastq . iterate ( options . stdin ) : #NEWLINE# #INDENT# if re . match ( options . grep_pattern , record . seq ) : #NEWLINE# #INDENT# options . stdout . write ( ( '%s\n' % record ) ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'reverse-complement' ) : #NEWLINE# #INDENT# for record in Fastq . iterate ( options . stdin ) : #NEWLINE# #INDENT# record . seq = Genomics . complement ( record . seq ) #NEWLINE# record . quals = record . quals [ : : ( - 1 ) ] #NEWLINE# options . stdout . write ( ( '%s\n' % record ) ) #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'sample' ) : #NEWLINE# #INDENT# sample_threshold = min ( 1.0 , options . sample_size ) #NEWLINE# random . seed ( options . seed ) #NEWLINE# if options . pair : #NEWLINE# #INDENT# if ( not options . output_filename_pattern ) : #NEWLINE# #INDENT# raise ValueError ( 'please specify output filename pattern for second pair (--output-filename-pattern)' ) #NEWLINE# #UNINDENT# outfile1 = options . stdout #NEWLINE# outfile2 = IOTools . openFile ( options . output_filename_pattern , 'w' ) #NEWLINE# for ( record1 , record2 ) in itertools . izip ( Fastq . iterate ( options . stdin ) , Fastq . iterate ( IOTools . openFile ( options . pair ) ) ) : #NEWLINE# #INDENT# c . input += 1 #NEWLINE# if ( random . random ( ) <= sample_threshold ) : #NEWLINE# #INDENT# c . output += 1 #NEWLINE# outfile1 . write ( ( '%s\n' % record1 ) ) #NEWLINE# outfile2 . write ( ( '%s\n' % record2 ) ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# for record in Fastq . iterate ( options . stdin ) : #NEWLINE# #INDENT# c . input += 1 #NEWLINE# if ( random . random ( ) <= sample_threshold ) : #NEWLINE# #INDENT# c . output += 1 #NEWLINE# options . stdout . write ( ( '%s\n' % record ) ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'apply' ) : #NEWLINE# #INDENT# ids = set ( IOTools . readList ( IOTools . openFile ( options . apply ) ) ) #NEWLINE# for record in Fastq . iterate ( options . stdin ) : #NEWLINE# #INDENT# c . input += 1 #NEWLINE# if ( re . sub ( ' .*' , '' , record . identifier ) . strip ( ) in ids ) : #NEWLINE# #INDENT# c . output += 1 #NEWLINE# options . stdout . write ( ( '%s\n' % record ) ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'trim3' ) : #NEWLINE# #INDENT# trim3 = options . nbases #NEWLINE# for record in Fastq . iterate ( options . stdin ) : #NEWLINE# #INDENT# c . input += 1 #NEWLINE# record . trim ( trim3 ) #NEWLINE# options . stdout . write ( ( '%s\n' % record ) ) #NEWLINE# c . output += 1 #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'trim5' ) : #NEWLINE# #INDENT# trim5 = options . nbases #NEWLINE# for record in Fastq . iterate ( options . stdin ) : #NEWLINE# #INDENT# c . input += 1 #NEWLINE# record . trim5 ( trim5 ) #NEWLINE# options . stdout . write ( ( '%s\n' % record ) ) #NEWLINE# c . output += 1 #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'unique' ) : #NEWLINE# #INDENT# keys = set ( ) #NEWLINE# for record in Fastq . iterate ( options . stdin ) : #NEWLINE# #INDENT# c . input += 1 #NEWLINE# if ( record . identifier in keys ) : #NEWLINE# #INDENT# continue #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# keys . add ( record . identifier ) #NEWLINE# #UNINDENT# options . stdout . write ( ( '%s\n' % record ) ) #NEWLINE# c . output += 1 #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'sort' ) : #NEWLINE# #INDENT# if ( not options . pair ) : #NEWLINE# #INDENT# statement = "paste - - - - | sort -k1,1 -t ' ' | tr '\t' '\n'" #NEWLINE# os . system ( statement ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( not options . output_filename_pattern ) : #NEWLINE# #INDENT# raise ValueError ( 'please specify output filename for second pair (--output-filename-pattern)' ) #NEWLINE# #UNINDENT# E . warn ( 'consider sorting individual fastq files - this is memory intensive' ) #NEWLINE# entries1 = { #NEWLINE# #INDENT# #NEWLINE# #UNINDENT# } #NEWLINE# entries2 = { #NEWLINE# #INDENT# #NEWLINE# #UNINDENT# } #NEWLINE# for ( record1 , record2 ) in itertools . izip ( Fastq . iterate ( options . stdin ) , Fastq . iterate ( IOTools . openFile ( options . pair ) ) ) : #NEWLINE# #INDENT# entries1 [ record1 . identifier [ : ( - 2 ) ] ] = ( record1 . seq , record1 . quals ) #NEWLINE# entries2 [ record2 . identifier [ : ( - 2 ) ] ] = ( record2 . seq , record2 . quals ) #NEWLINE# #UNINDENT# outfile1 = options . stdout #NEWLINE# outfile2 = IOTools . openFile ( options . output_filename_pattern , 'w' ) #NEWLINE# assert ( len ( set ( entries1 . keys ( ) ) . intersection ( set ( entries2 . keys ( ) ) ) ) == len ( entries1 ) ) , 'paired files do not contain the same reads need to reconcile files' #NEWLINE# for entry in sorted ( entries1 ) : #NEWLINE# #INDENT# outfile1 . write ( ( '@%s/1\n%s\n+\n%s\n' % ( entry , entries1 [ entry ] [ 0 ] , entries1 [ entry ] [ 1 ] ) ) ) #NEWLINE# outfile2 . write ( ( '@%s/2\n%s\n+\n%s\n' % ( entry , entries2 [ entry ] [ 0 ] , entries2 [ entry ] [ 1 ] ) ) ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'renumber-reads' ) : #NEWLINE# #INDENT# id_count = 1 #NEWLINE# for record in Fastq . iterate ( options . stdin ) : #NEWLINE# #INDENT# record . identifier = ( options . renumber_pattern % id_count ) #NEWLINE# id_count += 1 #NEWLINE# options . stdout . write ( ( '@%s\n%s\n+\n%s\n' % ( record . identifier , record . seq , record . quals ) ) ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# E . info ( ( '%s' % str ( c ) ) ) #NEWLINE# E . Stop ( )