root/galaxy-central/tools/sr_mapping/srma_wrapper.py

リビジョン 2, 10.6 KB (コミッタ: hatakeyama, 14 年 前)

import galaxy-central
行番号 
1#!/usr/bin/env python
2
3"""
4Runs SRMA on a SAM/BAM file;
5TODO: more documentation
6TODO: Allow use of fastq file in user's history, will need to index on-the-fly
7
8usage: 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
25import optparse, os, shutil, subprocess, sys, tempfile
26
27def stop_err( msg ):
28    sys.stderr.write( '%s\n' % msg )
29    sys.exit()
30
31def __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
202if __name__=="__main__": __main__()
Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。