Mercurial > repos > ryanmorin > nextgen_variant_identification
comparison SNV/snvmix.py @ 0:74f5ea818cea
Uploaded
author | ryanmorin |
---|---|
date | Wed, 12 Oct 2011 19:50:38 -0400 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:74f5ea818cea |
---|---|
1 #!/usr/bin/env python | |
2 | |
3 """ | |
4 Runs the SNVMix2 binary on a bam input file with various options. | |
5 | |
6 usage: %prog [options] | |
7 -i, --input1=i: bam file | |
8 -o, --output1=o: Output SNVMix (raw) | |
9 -d, --dbkey=d: dbkey of user-supplied file | |
10 -x, --indexDir=x: index directory | |
11 -t, --type=t: analysis type (e.g. mb|m|b|M|Mb|MB|SNVMix1) | |
12 -q, --base=q: base qual threshold | |
13 -Q, --map=Q: map qual threshold | |
14 -l, --pos=l: position file | |
15 -f, --full=f: Full mode (output scores for every position) | |
16 -R, --keep_dups: Retain reads flagged as duplicates (not recommended!) | |
17 -c, --keep_chastity: Retain reads that failed the chastity filter | |
18 """ | |
19 | |
20 import os, shutil, subprocess, sys, tempfile | |
21 from galaxy import eggs | |
22 import pkg_resources; pkg_resources.require( "bx-python" ) | |
23 from bx.cookbook import doc_optparse | |
24 | |
25 | |
26 | |
27 def stop_err( msg ): | |
28 sys.stderr.write( '%s\n' % msg ) | |
29 sys.exit() | |
30 | |
31 def check_seq_file( dbkey, GALAXY_DATA_INDEX_DIR ): | |
32 seqFile = '%s/sam_fa_indices.loc' % GALAXY_DATA_INDEX_DIR | |
33 seqPath = '' | |
34 for line in open( seqFile ): | |
35 line = line.rstrip( '\r\n' ) | |
36 if line and not line.startswith( '#' ) and line.startswith( 'index' ): | |
37 fields = line.split( '\t' ) | |
38 if len( fields ) < 3: | |
39 continue | |
40 if fields[1] == dbkey: | |
41 seqPath = fields[2].strip() | |
42 break | |
43 return seqPath | |
44 | |
45 def __main__(): | |
46 #Parse Command Line | |
47 options, args = doc_optparse.parse( __doc__ ) | |
48 seqPath = check_seq_file( options.dbkey, options.indexDir ) | |
49 | |
50 #make temp dir | |
51 tmpDir = tempfile.mkdtemp() | |
52 | |
53 #prepare basic SNVMix2 command | |
54 cmd = 'SNVMix2 -p b -i %s -r %s -o %s -q %s -Q %s -t %s' | |
55 try: | |
56 # have to nest try-except in try-finally to handle 2.4 | |
57 try: | |
58 if not os.path.exists( "%s.fai" % seqPath ): | |
59 raise Exception, "No sequences are available for '%s', request them by reporting this error." % options.dbkey | |
60 cmd = cmd % ( options.input1, seqPath, options.output1, options.base, options.map, options.type) | |
61 | |
62 if options.pos != "none": | |
63 if os.path.isfile(options.pos): | |
64 cmd = cmd + ' -l ' + options.pos | |
65 if options.full == "yes": | |
66 cmd = cmd + ' -f ' | |
67 else: | |
68 raise Exception, "position file doesn't exist" | |
69 | |
70 if options.keep_chastity == "yes": | |
71 cmd = cmd + ' -c' | |
72 if options.keep_dups == "yes": | |
73 cmd = cmd + ' -R' | |
74 | |
75 #perform SNVMix2 command | |
76 tmp = tempfile.NamedTemporaryFile( dir=tmpDir ).name | |
77 tmp_stderr = open( tmp, 'wb' ) | |
78 | |
79 proc = subprocess.Popen( args=cmd, shell=True, cwd=tmpDir, stderr=tmp_stderr.fileno() ) | |
80 returncode = proc.wait() | |
81 tmp_stderr.close() | |
82 #did it succeed? | |
83 # get stderr, allowing for case where it's very large | |
84 tmp_stderr = open( tmp, 'rb' ) | |
85 stderr = '' | |
86 buffsize = 1048576 | |
87 try: | |
88 while True: | |
89 stderr += tmp_stderr.read( buffsize ) | |
90 if not stderr or len( stderr ) % buffsize != 0: | |
91 break | |
92 except OverflowError: | |
93 pass | |
94 tmp_stderr.close() | |
95 if returncode != 0: | |
96 raise Exception, stderr | |
97 except Exception, e: | |
98 stop_err( 'Error running SNVMix2 tool\n' + str( e ) ) | |
99 finally: | |
100 #clean up temp files | |
101 if os.path.exists( tmpDir ): | |
102 shutil.rmtree( tmpDir ) | |
103 # check that there are results in the output file | |
104 if os.path.getsize( options.output1 ) > 0: | |
105 sys.stdout.write( 'wrote SNVMix output' ) | |
106 else: | |
107 stop_err( 'The output file is empty. Your input file may have had no matches, or there may be an error with your input file or settings.' ) | |
108 | |
109 if __name__ == "__main__" : __main__() |