Mercurial > repos > jjohnson > shear

annotate shear_wrapper.py @ 9:03071ff2e5c9

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add cap3 requirement for crest
author Jim Johnson <jj@umn.edu>
date Tue, 22 Oct 2013 09:55:49 -0500
parents 8fac6c7d32fd
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
4
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
1 #!/usr/bin/env python
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
2
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
3 import optparse, os, shutil, subprocess, sys, tempfile, os.path
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
4
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
5 def stop_err( msg, code ):
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
6 sys.stderr.write( '%s\n' % msg )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
7 sys.exit(code)
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
8
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
9 def __main__():
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
10 #Parse Command Line
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
11 parser = optparse.OptionParser()
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
12 parser.add_option( '-j', '--jar_path', dest='jar_path', help='Path to SHEAR.jar' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
13 parser.add_option( '-C', '--command', dest='command', help='SHEAR command to run' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
14 parser.add_option( '-p', '--prefix', dest='prefix', help='Prefix for all generated files. Required.' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
15 parser.add_option( '-b', '--bam', dest='bam', help='BAM alignment file containing the input sequences to the assembly.' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
16 parser.add_option( '-B', '--bai', dest='bai', help='BAM index file (.bai)' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
17 parser.add_option( '-f', '--fasta_reference', dest='fasta', help='The reference sequence fasta file' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
18 parser.add_option( '-F', '--fasta_index', dest='fai', help='The .fai index file for the reference sequence fasta' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
19 parser.add_option( '-t', '--twobit', dest='twobit', help='The .2bit encoding of the reference sequence fasta generated by faToTwoBit' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
20 parser.add_option( '-i', '--bwa_index', dest='bwa_index', help='The bwa index of the reference sequence' )
6
0158f7356ffd Add new outputs and params to shear_sv.xml for version 2.3
Jim Johnson <jj@umn.edu>
parents: 4
diff changeset
21 parser.add_option( '-R', '--report', dest='report', help='The SHEAR output report' )
4
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
22 parser.add_option( '-s', '--sdi', dest='sdi', help='The SHEAR sdi input from the SHEAR sv command' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
23 parser.add_option( '-o', '--output', dest='output', help='The SHEAR output assembly fasta file' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
24 parser.add_option( '-D', '--svidx_dir', dest='svidx_dir', help='The SHEAR output assembly fasta file' )
6
0158f7356ffd Add new outputs and params to shear_sv.xml for version 2.3
Jim Johnson <jj@umn.edu>
parents: 4
diff changeset
25 parser.add_option( '-S', '--sv-only', dest='sv_only', action="store_true", help='SV Only prediction mode.' )
0158f7356ffd Add new outputs and params to shear_sv.xml for version 2.3
Jim Johnson <jj@umn.edu>
parents: 4
diff changeset
26 parser.add_option( '-r', '--region', dest='region', help='Region of the input alignment to analyze' )
4
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
27 (options, args) = parser.parse_args()
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
28
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
29 def make_ref(src, dest, copy=False):
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
30 if copy:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
31 shutil.copy( src, dest )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
32 else:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
33 os.symlink( src, dest )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
34
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
35 # output version # of tool
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
36 #try:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
37 #except:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
38 # sys.stdout.write( 'Could not determine BWA version\n' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
39 if not options.jar_path:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
40 stop_err('path to SHEAR.jar is required',1)
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
41 if options.command:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
42 command = options.command
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
43 else:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
44 stop_err('SHEAR command is required',1)
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
45 args = [ 'java','-jar' ]
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
46 args.append( options.jar_path )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
47 args.append( command )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
48 # Check required params for command
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
49 buffsize = 1048576
8
8fac6c7d32fd Copy reference fasta for shear sv to avoid GATK issue with symlinks
Jim Johnson <jj@umn.edu>
parents: 6
diff changeset
50 copy_ref = True
4
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
51 prefix = options.prefix if options.prefix and len(options.prefix) > 0 else 'ref'
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
52 if command in ['sv']:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
53 args.append( '-p' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
54 args.append( prefix )
6
0158f7356ffd Add new outputs and params to shear_sv.xml for version 2.3
Jim Johnson <jj@umn.edu>
parents: 4
diff changeset
55 if options.sv_only:
0158f7356ffd Add new outputs and params to shear_sv.xml for version 2.3
Jim Johnson <jj@umn.edu>
parents: 4
diff changeset
56 args.append('--sv-only')
0158f7356ffd Add new outputs and params to shear_sv.xml for version 2.3
Jim Johnson <jj@umn.edu>
parents: 4
diff changeset
57 if options.region:
0158f7356ffd Add new outputs and params to shear_sv.xml for version 2.3
Jim Johnson <jj@umn.edu>
parents: 4
diff changeset
58 args.append('-r')
0158f7356ffd Add new outputs and params to shear_sv.xml for version 2.3
Jim Johnson <jj@umn.edu>
parents: 4
diff changeset
59 args.append(options.region)
8
8fac6c7d32fd Copy reference fasta for shear sv to avoid GATK issue with symlinks
Jim Johnson <jj@umn.edu>
parents: 6
diff changeset
60 copy_ref = True # GATK has issue with symlinks for .fa and .fai
4
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
61 if options.svidx_dir and command in ['sv']:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
62 if not os.path.isdir(options.svidx_dir):
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
63 os.makedirs(options.svidx_dir)
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
64 ref_prefix = os.path.join(options.svidx_dir,prefix)
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
65 copy_ref = True
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
66 else:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
67 ref_prefix = os.path.join(os.getcwd(),prefix)
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
68
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
69 if command in ['sv','assemble']:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
70 ref_fasta = ref_prefix + '.fa'
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
71 ref_index = ref_fasta + '.fai'
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
72 if options.fasta:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
73 make_ref( options.fasta, ref_fasta, copy=copy_ref )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
74 else:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
75 stop_err('Reference fasta file is required',3)
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
76 # if reference fasta index is not supllied, generate it
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
77 if options.fai:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
78 make_ref( options.fai, ref_index, copy=copy_ref )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
79 elif os.path.exists(options.fasta + '.fai'):
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
80 make_ref( options.fasta + '.fai', ref_index, copy=copy_ref )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
81 else:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
82 # generate fasta fai index
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
83 cmd1 = 'samtools faidx %s' % (ref_fasta )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
84 try:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
85 tmp_name = tempfile.NamedTemporaryFile(prefix='fai_', suffix='.err').name
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
86 tmp_stderr = open( tmp_name, 'wb' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
87 proc = subprocess.Popen( args=cmd1, shell=True, stderr=tmp_stderr.fileno() )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
88 returncode = proc.wait()
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
89 tmp_stderr.close()
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
90 if returncode != 0:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
91 stderr = ''
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
92 # get stderr, allowing for case where it's very large
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
93 tmp_stderr = open( tmp, 'rb' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
94 try:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
95 while True:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
96 stderr += tmp_stderr.read( buffsize )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
97 if not stderr or len( stderr ) % buffsize != 0:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
98 break
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
99 except OverflowError:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
100 pass
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
101 raise Exception, stderr
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
102 except Exception, e:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
103 stop_err( 'Error indexing reference sequence fasta file. ' + str( e ),5 )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
104 args.append( '-f' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
105 args.append( ref_fasta )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
106
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
107 if command in ['sv']:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
108 if not options.bam:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
109 stop_err('BAM alignment file is required',2)
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
110
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
111 # if reference 2bit is not supplied, generate it
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
112 twobit = ref_prefix + '.2bit'
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
113 if options.twobit:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
114 make_ref( options.twobit, twobit, copy=copy_ref )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
115 else:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
116 # generate 2bit index
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
117 cmd1 = 'faToTwoBit %s %s' % (ref_fasta, twobit )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
118 try:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
119 tmp_name = tempfile.NamedTemporaryFile(prefix='twobit_', suffix='.err').name
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
120 tmp_stderr = open( tmp_name, 'wb' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
121 proc = subprocess.Popen( args=cmd1, shell=True, stderr=tmp_stderr.fileno() )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
122 returncode = proc.wait()
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
123 tmp_stderr.close()
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
124 if returncode != 0:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
125 stderr = ''
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
126 # get stderr, allowing for case where it's very large
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
127 tmp_stderr = open( tmp, 'rb' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
128 try:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
129 while True:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
130 stderr += tmp_stderr.read( buffsize )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
131 if not stderr or len( stderr ) % buffsize != 0:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
132 break
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
133 except OverflowError:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
134 pass
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
135 raise Exception, stderr
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
136 except Exception, e:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
137 stop_err( 'Error indexing reference sequence fasta file. ' + str( e ),6 )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
138 args.append( '-t' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
139 args.append( twobit )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
140
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
141 # if bwa index is not supplied, generate it
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
142 bwa_index = ref_fasta
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
143 if options.bwa_index:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
144 if copy_ref:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
145 shutil.copytree(os.path.dirname(bwa_index),options.svidx_dir)
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
146 bwa_index = options.bwa_index
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
147 else:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
148 ONE_GB = 2**30
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
149 if os.stat( ref_fasta ).st_size <= ONE_GB: #use 1 GB as cut off for memory vs. max of 2gb database size; this is somewhat arbitrary
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
150 index_algorithm = 'is'
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
151 else:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
152 index_algorithm = 'bwtsw'
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
153 cmd1 = 'bwa index -a %s %s' % ( index_algorithm, ref_fasta )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
154 try:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
155 tmp_name = tempfile.NamedTemporaryFile(prefix='bwa_index_', suffix='.err').name
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
156 tmp_stderr = open( tmp_name, 'wb' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
157 proc = subprocess.Popen( args=cmd1, shell=True, stderr=tmp_stderr.fileno() )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
158 returncode = proc.wait()
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
159 tmp_stderr.close()
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
160 if returncode != 0:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
161 stderr = ''
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
162 # get stderr, allowing for case where it's very large
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
163 tmp_stderr = open( tmp, 'rb' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
164 try:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
165 while True:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
166 stderr += tmp_stderr.read( buffsize )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
167 if not stderr or len( stderr ) % buffsize != 0:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
168 break
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
169 except OverflowError:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
170 pass
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
171 raise Exception, stderr
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
172 bwa_index = ref_fasta
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
173 except Exception, e:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
174 # clean up temp dirs
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
175 stop_err( 'Error creating bwa index. ' + str( e ),7 )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
176 args.append( '-i' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
177 args.append( bwa_index )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
178
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
179 bam_file = os.path.join(os.getcwd(),prefix + '.bam')
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
180 bai_file = bam_file + '.bai'
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
181 if options.bam:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
182 os.symlink( options.bam, bam_file )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
183 else:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
184 stop_err('BAM alignment file is required',2)
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
185 # if bam index is not supplied, generate it
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
186 if options.bai:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
187 os.symlink( options.bai, bai_file )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
188 elif os.path.exists(options.bam + '.bai'):
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
189 os.symlink( options.bam + '.bai', bai_file )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
190 else:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
191 # generate bam index
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
192 cmd1 = 'samtools index %s %s' % (bam_file, bai_file )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
193 try:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
194 tmp_name = tempfile.NamedTemporaryFile(prefix='bai_', suffix='.err').name
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
195 tmp_stderr = open( tmp_name, 'wb' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
196 proc = subprocess.Popen( args=cmd1, shell=True, stderr=tmp_stderr.fileno() )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
197 returncode = proc.wait()
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
198 tmp_stderr.close()
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
199 if returncode != 0:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
200 stderr = ''
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
201 # get stderr, allowing for case where it's very large
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
202 tmp_stderr = open( tmp, 'rb' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
203 try:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
204 while True:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
205 stderr += tmp_stderr.read( buffsize )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
206 if not stderr or len( stderr ) % buffsize != 0:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
207 break
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
208 except OverflowError:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
209 pass
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
210 raise Exception, stderr
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
211 except Exception, e:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
212 stop_err( 'Error indexing BAM alignment file. ' + str( e ),4)
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
213 args.append( '-b' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
214 args.append( bam_file )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
215
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
216 if command in ['assemble']:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
217 if not options.sdi or not options.output:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
218 stop_err('SHEAR .sdi file and output file are required for assemble',2)
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
219 args.append( '-s' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
220 args.append( options.sdi )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
221 args.append( '-o' )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
222 args.append( options.output )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
223
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
224 # Run SHEAR command
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
225 print >> sys.stdout, "%s" % ' '.join(args)
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
226 try:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
227 proc = subprocess.Popen( args=args, shell=False )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
228 returncode = proc.wait()
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
229 if returncode != 0:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
230 stop_err( 'Error running SHEAR ' + command, returncode )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
231 except Exception, e:
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
232 # clean up temp dirs
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
233 stop_err( 'Error running SHEAR %s %s' % (command,str(e)),9 )
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
234
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
235 if __name__=="__main__": __main__()
a82400332451 Use shear_wrapper.py to generate reference indexes when needed
Jim Johnson <jj@umn.edu>
parents:
diff changeset
236