annotate SMART/DiffExpAnal/bam_to_sam_parallel.py @ 31:0ab839023fe4

Uploaded
author m-zytnicki
date Tue, 30 Apr 2013 14:33:21 -0400
parents 94ab73e8a190
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
18
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
1 #!/usr/bin/env python
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
2 """
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
3 Converts BAM data to sorted SAM data.
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
4 usage: bam_to_sam.py [options]
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
5 --input1: SAM file to be converted
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
6 --output1: output dataset in bam format
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
7 """
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
8
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
9 import optparse, os, sys, subprocess, tempfile, shutil, tarfile, random
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
10 #from galaxy import eggs
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
11 #import pkg_resources; pkg_resources.require( "bx-python" )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
12 #from bx.cookbook import doc_optparse
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
13 #from galaxy import util
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
14
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
15 def stop_err( msg ):
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
16 sys.stderr.write( '%s\n' % msg )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
17 sys.exit()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
18
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
19 def toTar(tarFileName, samOutputNames):
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
20 dir = os.path.dirname(tarFileName)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
21 tfile = tarfile.open(tarFileName + ".tmp.tar", "w")
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
22 currentPath = os.getcwd()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
23 os.chdir(dir)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
24 for file in samOutputNames:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
25 relativeFileName = os.path.basename(file)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
26 tfile.add(relativeFileName)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
27 os.system("mv %s %s" % (tarFileName + ".tmp.tar", tarFileName))
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
28 tfile.close()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
29 os.chdir(currentPath)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
30
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
31
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
32 def __main__():
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
33 #Parse Command Line
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
34 parser = optparse.OptionParser()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
35 parser.add_option('-t', '--tar', dest='outputTar', default=None, help='output all SAM results in a tar file.' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
36 parser.add_option( '', '--input1', dest='input1', help='The input list of BAM datasets on txt format.' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
37 #parser.add_option( '', '--input1', dest='input1', help='The input BAM dataset' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
38 parser.add_option( '', '--output1', dest='output1', help='The output list of SAM datasets on txt format.' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
39 #parser.add_option( '', '--output1', dest='output1', help='The output SAM dataset' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
40 parser.add_option( '', '--header', dest='header', action='store_true', default=False, help='Write SAM Header' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
41 ( options, args ) = parser.parse_args()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
42
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
43
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
44 #Parse the input txt file and read a list of BAM files.
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
45 file = open(options.input1, "r")
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
46 lines = file.readlines()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
47 inputFileNames = []
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
48 samOutputNames = []
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
49 outputName = options.output1
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
50 resDirName = os.path.dirname(outputName) + '/'
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
51 #Write output txt file and define all output sam file names.
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
52 out = open(outputName, "w")
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
53 for line in lines:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
54 tab = line.split()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
55 inputFileNames.append(tab[1])
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
56 samOutName = resDirName + tab[0] + '_samOutput_%s.sam' % random.randrange(0, 10000)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
57 samOutputNames.append(samOutName)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
58 out.write(tab[0] + '\t' + samOutName + '\n')
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
59 file.close()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
60 out.close()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
61
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
62 # output version # of tool
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
63 try:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
64 tmp_files = []
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
65 tmp = tempfile.NamedTemporaryFile().name
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
66 tmp_files.append(tmp)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
67 tmp_stdout = open( tmp, 'wb' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
68 proc = subprocess.Popen( args='samtools 2>&1', shell=True, stdout=tmp_stdout )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
69 tmp_stdout.close()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
70 returncode = proc.wait()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
71 stdout = None
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
72 for line in open( tmp_stdout.name, 'rb' ):
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
73 if line.lower().find( 'version' ) >= 0:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
74 stdout = line.strip()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
75 break
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
76 if stdout:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
77 sys.stdout.write( 'Samtools %s\n' % stdout )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
78 else:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
79 raise Exception
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
80 except:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
81 sys.stdout.write( 'Could not determine Samtools version\n' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
82
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
83
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
84
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
85 tmp_dirs = []
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
86 for i in range(len(inputFileNames)):
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
87 try:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
88 # exit if input file empty
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
89 if os.path.getsize( inputFileNames[i] ) == 0:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
90 raise Exception, 'Initial input txt file is empty.'
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
91 # Sort alignments by leftmost coordinates. File <out.prefix>.bam will be created. This command
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
92 # may also create temporary files <out.prefix>.%d.bam when the whole alignment cannot be fitted
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
93 # into memory ( controlled by option -m ).
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
94 tmp_dir = tempfile.mkdtemp()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
95 tmp_dirs.append(tmp_dir)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
96 tmp_sorted_aligns_file = tempfile.NamedTemporaryFile( dir=tmp_dir )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
97 tmp_sorted_aligns_file_base = tmp_sorted_aligns_file.name
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
98 tmp_sorted_aligns_file_name = '%s.bam' % tmp_sorted_aligns_file.name
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
99 tmp_files.append(tmp_sorted_aligns_file_name)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
100 tmp_sorted_aligns_file.close()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
101
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
102 command = 'samtools sort %s %s' % ( inputFileNames[i], tmp_sorted_aligns_file_base )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
103 tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
104 tmp_stderr = open( tmp, 'wb' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
105 proc = subprocess.Popen( args=command, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
106 returncode = proc.wait()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
107 tmp_stderr.close()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
108 # get stderr, allowing for case where it's very large
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
109 tmp_stderr = open( tmp, 'rb' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
110 stderr = ''
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
111 buffsize = 1048576
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
112 try:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
113 while True:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
114 stderr += tmp_stderr.read( buffsize )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
115 if not stderr or len( stderr ) % buffsize != 0:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
116 break
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
117 except OverflowError:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
118 pass
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
119 tmp_stderr.close()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
120 if returncode != 0:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
121 raise Exception, stderr
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
122 # exit if sorted BAM file empty
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
123 if os.path.getsize( tmp_sorted_aligns_file_name) == 0:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
124 raise Exception, 'Intermediate sorted BAM file empty'
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
125 except Exception, e:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
126 stop_err( 'Error sorting alignments from (%s), %s' % ( inputFileNames[i], str( e ) ) )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
127
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
128 try:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
129 # Extract all alignments from the input BAM file to SAM format ( since no region is specified, all the alignments will be extracted ).
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
130 if options.header:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
131 view_options = "-h"
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
132 else:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
133 view_options = ""
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
134 command = 'samtools view %s -o %s %s' % ( view_options, samOutputNames[i], tmp_sorted_aligns_file_name )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
135 tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
136 tmp_stderr = open( tmp, 'wb' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
137 proc = subprocess.Popen( args=command, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
138 returncode = proc.wait()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
139 tmp_stderr.close()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
140 # get stderr, allowing for case where it's very large
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
141 tmp_stderr = open( tmp, 'rb' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
142 stderr = ''
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
143 buffsize = 1048576
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
144 try:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
145 while True:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
146 stderr += tmp_stderr.read( buffsize )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
147 if not stderr or len( stderr ) % buffsize != 0:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
148 break
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
149 except OverflowError:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
150 pass
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
151 tmp_stderr.close()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
152 if returncode != 0:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
153 raise Exception, stderr
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
154 except Exception, e:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
155 stop_err( 'Error extracting alignments from (%s), %s' % ( inputFileNames[i], str( e ) ) )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
156 if os.path.getsize( samOutputNames[i] ) > 0:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
157 sys.stdout.write( 'BAM file converted to SAM' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
158 else:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
159 stop_err( 'The output file is empty, there may be an error with your input file.' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
160
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
161 if options.outputTar != None:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
162 toTar(options.outputTar, samOutputNames)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
163 #clean up temp files
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
164 for tmp_dir in tmp_dirs:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
165 if os.path.exists( tmp_dir ):
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
166 shutil.rmtree( tmp_dir )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
167 #print tmp_files
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
168 #for tmp in tmp_files:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
169 # os.remove(tmp)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
170
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
171
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
172 if __name__=="__main__": __main__()