annotate bam_to_scidx.py @ 1:24a8165055bd draft default tip

Uploaded
author greg
date Wed, 02 Dec 2015 16:13:17 -0500
parents 962b0366a980
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
962b0366a980 Uploaded
greg
parents:
diff changeset
1 """
962b0366a980 Uploaded
greg
parents:
diff changeset
2 bam_to_scidx.py
962b0366a980 Uploaded
greg
parents:
diff changeset
3
962b0366a980 Uploaded
greg
parents:
diff changeset
4 Input: BAM file
962b0366a980 Uploaded
greg
parents:
diff changeset
5 Output: Converted scidx file
962b0366a980 Uploaded
greg
parents:
diff changeset
6 """
962b0366a980 Uploaded
greg
parents:
diff changeset
7 import optparse
962b0366a980 Uploaded
greg
parents:
diff changeset
8 import os
962b0366a980 Uploaded
greg
parents:
diff changeset
9 import shutil
962b0366a980 Uploaded
greg
parents:
diff changeset
10 import subprocess
962b0366a980 Uploaded
greg
parents:
diff changeset
11 import sys
962b0366a980 Uploaded
greg
parents:
diff changeset
12 import tempfile
962b0366a980 Uploaded
greg
parents:
diff changeset
13
962b0366a980 Uploaded
greg
parents:
diff changeset
14 BUFF_SIZE = 1048576
962b0366a980 Uploaded
greg
parents:
diff changeset
15
962b0366a980 Uploaded
greg
parents:
diff changeset
16
962b0366a980 Uploaded
greg
parents:
diff changeset
17 def get_stderr_exception(tmp_err, tmp_stderr):
962b0366a980 Uploaded
greg
parents:
diff changeset
18 """
962b0366a980 Uploaded
greg
parents:
diff changeset
19 Return a stderr string of reasonable size.
962b0366a980 Uploaded
greg
parents:
diff changeset
20 """
962b0366a980 Uploaded
greg
parents:
diff changeset
21 tmp_stderr.close()
962b0366a980 Uploaded
greg
parents:
diff changeset
22 # Get stderr, allowing for case where it's very large.
962b0366a980 Uploaded
greg
parents:
diff changeset
23 tmp_stderr = open(tmp_err, 'rb')
962b0366a980 Uploaded
greg
parents:
diff changeset
24 stderr_str = ''
962b0366a980 Uploaded
greg
parents:
diff changeset
25 buffsize = BUFF_SIZE
962b0366a980 Uploaded
greg
parents:
diff changeset
26 try:
962b0366a980 Uploaded
greg
parents:
diff changeset
27 while True:
962b0366a980 Uploaded
greg
parents:
diff changeset
28 stderr_str += tmp_stderr.read(buffsize)
962b0366a980 Uploaded
greg
parents:
diff changeset
29 if not stderr_str or len(stderr_str) % buffsize != 0:
962b0366a980 Uploaded
greg
parents:
diff changeset
30 break
962b0366a980 Uploaded
greg
parents:
diff changeset
31 except OverflowError:
962b0366a980 Uploaded
greg
parents:
diff changeset
32 pass
962b0366a980 Uploaded
greg
parents:
diff changeset
33 tmp_stderr.close()
962b0366a980 Uploaded
greg
parents:
diff changeset
34 return stderr_str
962b0366a980 Uploaded
greg
parents:
diff changeset
35
962b0366a980 Uploaded
greg
parents:
diff changeset
36
962b0366a980 Uploaded
greg
parents:
diff changeset
37 def stop_err(msg):
962b0366a980 Uploaded
greg
parents:
diff changeset
38 sys.stderr.write(msg)
962b0366a980 Uploaded
greg
parents:
diff changeset
39 sys.exit(1)
962b0366a980 Uploaded
greg
parents:
diff changeset
40
962b0366a980 Uploaded
greg
parents:
diff changeset
41
962b0366a980 Uploaded
greg
parents:
diff changeset
42 parser = optparse.OptionParser()
962b0366a980 Uploaded
greg
parents:
diff changeset
43 parser.add_option('-j', '--jar_file', dest='jar_file', type='string', help='BAMtoIDX.jar')
962b0366a980 Uploaded
greg
parents:
diff changeset
44 parser.add_option('-b', '--input_bam', dest='input_bam', type='string', help='Input dataset in BAM format')
962b0366a980 Uploaded
greg
parents:
diff changeset
45 parser.add_option('-i', '--input_bai', dest='input_bai', type='string', help='Input dataset index')
962b0366a980 Uploaded
greg
parents:
diff changeset
46 parser.add_option('-r', '--read', dest='read', type='int', default=0, help='Reads.')
962b0366a980 Uploaded
greg
parents:
diff changeset
47 parser.add_option('-o', '--output', dest='output', type='string', help='Output dataset in scidx format')
962b0366a980 Uploaded
greg
parents:
diff changeset
48 options, args = parser.parse_args()
962b0366a980 Uploaded
greg
parents:
diff changeset
49
962b0366a980 Uploaded
greg
parents:
diff changeset
50 tmp_dir = tempfile.mkdtemp(prefix='tmp-scidx-')
962b0366a980 Uploaded
greg
parents:
diff changeset
51 tmp_out = tempfile.NamedTemporaryFile().name
962b0366a980 Uploaded
greg
parents:
diff changeset
52 tmp_stdout = open(tmp_out, 'wb')
962b0366a980 Uploaded
greg
parents:
diff changeset
53 tmp_err = tempfile.NamedTemporaryFile().name
962b0366a980 Uploaded
greg
parents:
diff changeset
54 tmp_stderr = open(tmp_err, 'wb')
962b0366a980 Uploaded
greg
parents:
diff changeset
55
962b0366a980 Uploaded
greg
parents:
diff changeset
56 # Link input BAM file and associated index file into the working directory.
962b0366a980 Uploaded
greg
parents:
diff changeset
57 input_data_file_name = "input.bam"
962b0366a980 Uploaded
greg
parents:
diff changeset
58 os.link(options.input_bam, os.path.join(tmp_dir, input_data_file_name))
962b0366a980 Uploaded
greg
parents:
diff changeset
59 os.link(options.input_bai, os.path.join(tmp_dir, "%s.bai" % input_data_file_name))
962b0366a980 Uploaded
greg
parents:
diff changeset
60
962b0366a980 Uploaded
greg
parents:
diff changeset
61 cmd = "java -jar %s -b %s -s %d -o %s" % (options.jar_file, input_data_file_name, options.read, options.output)
962b0366a980 Uploaded
greg
parents:
diff changeset
62 proc = subprocess.Popen(args=cmd, stdout=tmp_stdout, stderr=tmp_stderr, shell=True, cwd=tmp_dir)
962b0366a980 Uploaded
greg
parents:
diff changeset
63 return_code = proc.wait()
962b0366a980 Uploaded
greg
parents:
diff changeset
64 if return_code != 0:
962b0366a980 Uploaded
greg
parents:
diff changeset
65 error_message = get_stderr_exception(tmp_err, tmp_stderr)
962b0366a980 Uploaded
greg
parents:
diff changeset
66 stop_err(error_message)
962b0366a980 Uploaded
greg
parents:
diff changeset
67
962b0366a980 Uploaded
greg
parents:
diff changeset
68 if tmp_dir and os.path.exists(tmp_dir):
962b0366a980 Uploaded
greg
parents:
diff changeset
69 shutil.rmtree(tmp_dir)