[2] | 1 | #!/usr/bin/env python |
---|
| 2 | |
---|
| 3 | """ |
---|
| 4 | Runs SRMA on a SAM/BAM file; |
---|
| 5 | TODO: more documentation |
---|
| 6 | TODO: Allow use of fastq file in user's history, will need to index on-the-fly |
---|
| 7 | |
---|
| 8 | usage: srma_wrapper.py [options] |
---|
| 9 | -r, --ref=r: The reference genome to use or index |
---|
| 10 | -i, --input=i: The SAM/BAM input file |
---|
| 11 | -o, --output=o: The SAM/BAM output file |
---|
| 12 | -O, --offset=O: The alignment offset |
---|
| 13 | -Q, --minMappingQuality=Q: The minimum mapping quality |
---|
| 14 | -P, --minAlleleProbability=P: The minimum allele probability conditioned on coverage (for the binomial quantile). |
---|
| 15 | -C, --minAlleleCoverage=C: The minimum haploid coverage for the consensus. Default value: 3. This option can be set |
---|
| 16 | -R, --range=R: A range to examine |
---|
| 17 | -c, --correctBases=c: Correct bases |
---|
| 18 | -q, --useSequenceQualities=q: Use sequence qualities |
---|
| 19 | -M, --maxHeapSize=M: The maximum number of nodes on the heap before re-alignment is ignored |
---|
| 20 | -s, --fileSource=s: Whether to use a previously indexed reference sequence or one from history (indexed or history) |
---|
| 21 | -p, --params=p: Parameter setting to use (pre_set or full) |
---|
| 22 | -D, --dbkey=D: Dbkey for reference genome |
---|
| 23 | """ |
---|
| 24 | |
---|
| 25 | import optparse, os, shutil, subprocess, sys, tempfile |
---|
| 26 | |
---|
| 27 | def stop_err( msg ): |
---|
| 28 | sys.stderr.write( '%s\n' % msg ) |
---|
| 29 | sys.exit() |
---|
| 30 | |
---|
| 31 | def __main__(): |
---|
| 32 | #Parse Command Line |
---|
| 33 | parser = optparse.OptionParser() |
---|
| 34 | parser.add_option( '-r', '--ref', dest='ref', help='The reference genome to index and use' ) |
---|
| 35 | parser.add_option( '-u', '--refUID', dest='refUID', help='The pre-index reference genome unique Identifier' ) |
---|
| 36 | parser.add_option( '-L', '--refLocations', dest='refLocations', help='The filepath to the srma indices location file' ) |
---|
| 37 | parser.add_option( '-i', '--input', dest='input', help='The SAM/BAM input file' ) |
---|
| 38 | parser.add_option( '-I', '--inputIndex', dest='inputIndex', help='The SAM/BAM input index file' ) |
---|
| 39 | parser.add_option( '-o', '--output', dest='output', help='The SAM/BAM output file' ) |
---|
| 40 | parser.add_option( '-O', '--offset', dest='offset', help='The alignment offset' ) |
---|
| 41 | parser.add_option( '-Q', '--minMappingQuality', dest='minMappingQuality', help='The minimum mapping quality' ) |
---|
| 42 | parser.add_option( '-P', '--minAlleleProbability', dest='minAlleleProbability', help='The minimum allele probability conditioned on coverage (for the binomial quantile).' ) |
---|
| 43 | parser.add_option( '-C', '--minAlleleCoverage', dest='minAlleleCoverage', help='The minimum haploid coverage for the consensus' ) |
---|
| 44 | parser.add_option( '-R', '--range', dest='range', help='A range to examine' ) |
---|
| 45 | parser.add_option( '-c', '--correctBases', dest='correctBases', help='Correct bases ' ) |
---|
| 46 | parser.add_option( '-q', '--useSequenceQualities', dest='useSequenceQualities', help='Use sequence qualities ' ) |
---|
| 47 | parser.add_option( '-M', '--maxHeapSize', dest='maxHeapSize', help='The maximum number of nodes on the heap before re-alignment is ignored' ) |
---|
| 48 | parser.add_option( '-s', '--fileSource', dest='fileSource', help='Whether to use a previously indexed reference sequence or one from history (indexed or history)' ) |
---|
| 49 | parser.add_option( '-p', '--params', dest='params', help='Parameter setting to use (pre_set or full)' ) |
---|
| 50 | parser.add_option( '-j', '--jarBin', dest='jarBin', default='', help='The path to where jars are stored' ) |
---|
| 51 | parser.add_option( '-f', '--jarFile', dest='jarFile', help='The file name of the jar file to use') |
---|
| 52 | (options, args) = parser.parse_args() |
---|
| 53 | |
---|
| 54 | # make temp directory for srma |
---|
| 55 | tmp_dir = tempfile.mkdtemp() |
---|
| 56 | buffsize = 1048576 |
---|
| 57 | |
---|
| 58 | # set up reference filenames |
---|
| 59 | reference_filepath_name = None |
---|
| 60 | # need to create SRMA dict and Samtools fai files for custom genome |
---|
| 61 | if options.fileSource == 'history': |
---|
| 62 | try: |
---|
| 63 | reference_filepath = tempfile.NamedTemporaryFile( dir=tmp_dir, suffix='.fa' ) |
---|
| 64 | reference_filepath_name = reference_filepath.name |
---|
| 65 | reference_filepath.close() |
---|
| 66 | fai_filepath_name = '%s.fai' % reference_filepath_name |
---|
| 67 | dict_filepath_name = reference_filepath_name.replace( '.fa', '.dict' ) |
---|
| 68 | os.symlink( options.ref, reference_filepath_name ) |
---|
| 69 | # create fai file using Samtools |
---|
| 70 | index_fai_cmd = 'samtools faidx %s' % reference_filepath_name |
---|
| 71 | try: |
---|
| 72 | tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name |
---|
| 73 | tmp_stderr = open( tmp, 'wb' ) |
---|
| 74 | proc = subprocess.Popen( args=index_fai_cmd, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() ) |
---|
| 75 | returncode = proc.wait() |
---|
| 76 | tmp_stderr.close() |
---|
| 77 | # get stderr, allowing for case where it's very large |
---|
| 78 | tmp_stderr = open( tmp, 'rb' ) |
---|
| 79 | stderr = '' |
---|
| 80 | try: |
---|
| 81 | while True: |
---|
| 82 | stderr += tmp_stderr.read( buffsize ) |
---|
| 83 | if not stderr or len( stderr ) % buffsize != 0: |
---|
| 84 | break |
---|
| 85 | except OverflowError: |
---|
| 86 | pass |
---|
| 87 | tmp_stderr.close() |
---|
| 88 | if returncode != 0: |
---|
| 89 | raise Exception, stderr |
---|
| 90 | except Exception, e: |
---|
| 91 | # clean up temp dir |
---|
| 92 | if os.path.exists( tmp_dir ): |
---|
| 93 | shutil.rmtree( tmp_dir ) |
---|
| 94 | stop_err( 'Error creating Samtools index for custom genome file: %s\n' % str( e ) ) |
---|
| 95 | # create dict file using SRMA |
---|
| 96 | dict_cmd = 'java -cp "%s" net.sf.picard.sam.CreateSequenceDictionary R=%s O=%s' % ( os.path.join( options.jarBin, options.jarFile ), reference_filepath_name, dict_filepath_name ) |
---|
| 97 | try: |
---|
| 98 | tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name |
---|
| 99 | tmp_stderr = open( tmp, 'wb' ) |
---|
| 100 | proc = subprocess.Popen( args=dict_cmd, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() ) |
---|
| 101 | returncode = proc.wait() |
---|
| 102 | tmp_stderr.close() |
---|
| 103 | # get stderr, allowing for case where it's very large |
---|
| 104 | tmp_stderr = open( tmp, 'rb' ) |
---|
| 105 | stderr = '' |
---|
| 106 | try: |
---|
| 107 | while True: |
---|
| 108 | stderr += tmp_stderr.read( buffsize ) |
---|
| 109 | if not stderr or len( stderr ) % buffsize != 0: |
---|
| 110 | break |
---|
| 111 | except OverflowError: |
---|
| 112 | pass |
---|
| 113 | tmp_stderr.close() |
---|
| 114 | if returncode != 0: |
---|
| 115 | raise Exception, stderr |
---|
| 116 | except Exception, e: |
---|
| 117 | # clean up temp dir |
---|
| 118 | if os.path.exists( tmp_dir ): |
---|
| 119 | shutil.rmtree( tmp_dir ) |
---|
| 120 | stop_err( 'Error creating index for custom genome file: %s\n' % str( e ) ) |
---|
| 121 | except Exception, e: |
---|
| 122 | # clean up temp dir |
---|
| 123 | if os.path.exists( tmp_dir ): |
---|
| 124 | shutil.rmtree( tmp_dir ) |
---|
| 125 | stop_err( 'Problem handling SRMA index (dict file) for custom genome file: %s\n' % str( e ) ) |
---|
| 126 | # using built-in dict/index files |
---|
| 127 | else: |
---|
| 128 | for line in open( options.refLocations ): |
---|
| 129 | if not line.startswith( '#' ): |
---|
| 130 | fields = line.strip().split( '\t' ) |
---|
| 131 | if len( fields ) >= 2: |
---|
| 132 | if fields[0] == options.refUID: |
---|
| 133 | reference_filepath_name = fields[1] |
---|
| 134 | break |
---|
| 135 | if reference_filepath_name is None: |
---|
| 136 | raise ValueError( 'A valid genome reference was not provided.' ) |
---|
| 137 | |
---|
| 138 | # set up aligning and generate aligning command options |
---|
| 139 | if options.params == 'pre_set': |
---|
| 140 | srma_cmds = '' |
---|
| 141 | else: |
---|
| 142 | if options.useSequenceQualities == 'true': |
---|
| 143 | useSequenceQualities = 'true' |
---|
| 144 | else: |
---|
| 145 | useSequenceQualities = 'false' |
---|
| 146 | ranges = 'null' |
---|
| 147 | if options.range == 'None': |
---|
| 148 | range = 'null' |
---|
| 149 | else: |
---|
| 150 | range = options.range |
---|
| 151 | srma_cmds = "OFFSET=%s MIN_MAPQ=%s MINIMUM_ALLELE_PROBABILITY=%s MINIMUM_ALLELE_COVERAGE=%s RANGES=%s RANGE=%s CORRECT_BASES=%s USE_SEQUENCE_QUALITIES=%s MAX_HEAP_SIZE=%s" % ( options.offset, options.minMappingQuality, options.minAlleleProbability, options.minAlleleCoverage, ranges, range, options.correctBases, options.useSequenceQualities, options.maxHeapSize ) |
---|
| 152 | |
---|
| 153 | # perform alignments |
---|
| 154 | buffsize = 1048576 |
---|
| 155 | try: |
---|
| 156 | #symlink input bam and index files due to the naming conventions required by srma here |
---|
| 157 | input_bam_filename = os.path.join( tmp_dir, '%s.bam' % os.path.split( options.input )[-1] ) |
---|
| 158 | os.symlink( options.input, input_bam_filename ) |
---|
| 159 | input_bai_filename = "%s.bai" % os.path.splitext( input_bam_filename )[0] |
---|
| 160 | os.symlink( options.inputIndex, input_bai_filename ) |
---|
| 161 | |
---|
| 162 | #create a temp output name, ending in .bam due to required naming conventions? unkown if required |
---|
| 163 | output_bam_filename = os.path.join( tmp_dir, "%s.bam" % os.path.split( options.output )[-1] ) |
---|
| 164 | # generate commandline |
---|
| 165 | cmd = 'java -jar %s I=%s O=%s R=%s %s' % ( os.path.join( options.jarBin, options.jarFile ), input_bam_filename, output_bam_filename, reference_filepath_name, srma_cmds ) |
---|
| 166 | |
---|
| 167 | # need to nest try-except in try-finally to handle 2.4 |
---|
| 168 | try: |
---|
| 169 | try: |
---|
| 170 | tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name |
---|
| 171 | tmp_stderr = open( tmp, 'wb' ) |
---|
| 172 | proc = subprocess.Popen( args=cmd, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() ) |
---|
| 173 | returncode = proc.wait() |
---|
| 174 | tmp_stderr.close() |
---|
| 175 | # get stderr, allowing for case where it's very large |
---|
| 176 | tmp_stderr = open( tmp, 'rb' ) |
---|
| 177 | stderr = '' |
---|
| 178 | try: |
---|
| 179 | while True: |
---|
| 180 | stderr += tmp_stderr.read( buffsize ) |
---|
| 181 | if not stderr or len( stderr ) % buffsize != 0: |
---|
| 182 | break |
---|
| 183 | except OverflowError: |
---|
| 184 | pass |
---|
| 185 | tmp_stderr.close() |
---|
| 186 | if returncode != 0: |
---|
| 187 | raise Exception, stderr |
---|
| 188 | except Exception, e: |
---|
| 189 | raise Exception, 'Error executing SRMA. ' + str( e ) |
---|
| 190 | # move file from temp location (with .bam name) to provided path |
---|
| 191 | shutil.move( output_bam_filename, options.output ) |
---|
| 192 | # check that there are results in the output file |
---|
| 193 | if os.path.getsize( options.output ) <= 0: |
---|
| 194 | raise Exception, 'The output file is empty. You may simply have no matches, or there may be an error with your input file or settings.' |
---|
| 195 | except Exception, e: |
---|
| 196 | stop_err( 'The re-alignment failed.\n' + str( e ) ) |
---|
| 197 | finally: |
---|
| 198 | # clean up temp dir |
---|
| 199 | if os.path.exists( tmp_dir ): |
---|
| 200 | shutil.rmtree( tmp_dir ) |
---|
| 201 | |
---|
| 202 | if __name__=="__main__": __main__() |
---|