Mercurial > repos > ryanmorin > nextgen_variant_identification

comparison SNV/snvmix.py @ 0:74f5ea818cea

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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__()