annotate sam_to_bam.py @ 5:c73bf16b45df draft

Uploaded
author devteam
date Thu, 05 Mar 2015 21:28:25 -0500
parents ab4c4e07eb3c
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
1 #!/usr/bin/env python
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
2 """
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
3 Converts SAM data to sorted BAM data.
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
4 usage: sam_to_bam.py [options]
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
5 --input1: SAM file to be converted
1
93f2e3337a33 Update sam_to_bam to use the fasta_indexes data table.
Dave Bouvier <dave@bx.psu.edu>
parents: 0
diff changeset
6 --index: path of the indexed reference genome
0
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
7 --ref_file: Reference file if choosing from history
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
8 --output1: output dataset in bam format
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
9 """
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
10
1
93f2e3337a33 Update sam_to_bam to use the fasta_indexes data table.
Dave Bouvier <dave@bx.psu.edu>
parents: 0
diff changeset
11 import optparse, os, sys, subprocess, tempfile, shutil
0
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
12
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
13 def stop_err( msg ):
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
14 sys.stderr.write( '%s\n' % msg )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
15 sys.exit()
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
16
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
17 def __main__():
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
18 #Parse Command Line
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
19 parser = optparse.OptionParser()
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
20 parser.add_option( '', '--input1', dest='input1', help='The input SAM dataset' )
1
93f2e3337a33 Update sam_to_bam to use the fasta_indexes data table.
Dave Bouvier <dave@bx.psu.edu>
parents: 0
diff changeset
21
93f2e3337a33 Update sam_to_bam to use the fasta_indexes data table.
Dave Bouvier <dave@bx.psu.edu>
parents: 0
diff changeset
22 parser.add_option( '', '--index', dest='index', help='The path of the indexed reference genome' )
0
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
23 parser.add_option( '', '--ref_file', dest='ref_file', help='The reference dataset from the history' )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
24 parser.add_option( '', '--output1', dest='output1', help='The output BAM dataset' )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
25 ( options, args ) = parser.parse_args()
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
26
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
27 # output version # of tool
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
28 try:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
29 tmp = tempfile.NamedTemporaryFile().name
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
30 tmp_stdout = open( tmp, 'wb' )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
31 proc = subprocess.Popen( args='samtools 2>&1', shell=True, stdout=tmp_stdout )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
32 tmp_stdout.close()
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
33 returncode = proc.wait()
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
34 stdout = None
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
35 for line in open( tmp_stdout.name, 'rb' ):
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
36 if line.lower().find( 'version' ) >= 0:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
37 stdout = line.strip()
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
38 break
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
39 if stdout:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
40 sys.stdout.write( 'Samtools %s\n' % stdout )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
41 else:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
42 raise Exception
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
43 except:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
44 sys.stdout.write( 'Could not determine Samtools version\n' )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
45
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
46 tmp_dir = tempfile.mkdtemp( dir='.' )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
47 if not options.ref_file or options.ref_file == 'None':
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
48 # We're using locally cached reference sequences( e.g., /galaxy/data/equCab2/sam_index/equCab2.fa ).
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
49 # The indexes for /galaxy/data/equCab2/sam_index/equCab2.fa will be contained in
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
50 # a file named /galaxy/data/equCab2/sam_index/equCab2.fa.fai
1
93f2e3337a33 Update sam_to_bam to use the fasta_indexes data table.
Dave Bouvier <dave@bx.psu.edu>
parents: 0
diff changeset
51 fai_index_file_path = '%s.fai' % options.index
0
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
52 if not os.path.exists( fai_index_file_path ):
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
53 #clean up temp files
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
54 if os.path.exists( tmp_dir ):
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
55 shutil.rmtree( tmp_dir )
1
93f2e3337a33 Update sam_to_bam to use the fasta_indexes data table.
Dave Bouvier <dave@bx.psu.edu>
parents: 0
diff changeset
56 stop_err( 'Indexed genome %s not present, request it by reporting this error.' % options.index )
0
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
57 else:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
58 try:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
59 # Create indexes for history reference ( e.g., ~/database/files/000/dataset_1.dat ) using samtools faidx, which will:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
60 # - index reference sequence in the FASTA format or extract subsequence from indexed reference sequence
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
61 # - if no region is specified, faidx will index the file and create <ref.fasta>.fai on the disk
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
62 # - if regions are specified, the subsequences will be retrieved and printed to stdout in the FASTA format
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
63 # - the input file can be compressed in the RAZF format.
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
64 # IMPORTANT NOTE: a real weakness here is that we are creating indexes for the history dataset
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
65 # every time we run this tool. It would be nice if we could somehow keep track of user's specific
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
66 # index files so they could be re-used.
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
67 fai_index_file_base = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
68 # At this point, fai_index_file_path will look something like /tmp/dataset_13.dat
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
69 os.symlink( options.ref_file, fai_index_file_base )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
70 fai_index_file_path = '%s.fai' % fai_index_file_base
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
71 command = 'samtools faidx %s' % fai_index_file_base
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
72 tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
73 tmp_stderr = open( tmp, 'wb' )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
74 proc = subprocess.Popen( args=command, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
75 returncode = proc.wait()
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
76 tmp_stderr.close()
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
77 # get stderr, allowing for case where it's very large
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
78 tmp_stderr = open( tmp, 'rb' )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
79 stderr = ''
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
80 buffsize = 1048576
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
81 try:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
82 while True:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
83 stderr += tmp_stderr.read( buffsize )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
84 if not stderr or len( stderr ) % buffsize != 0:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
85 break
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
86 except OverflowError:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
87 pass
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
88 tmp_stderr.close()
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
89 if returncode != 0:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
90 raise Exception, stderr
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
91 if os.path.getsize( fai_index_file_path ) == 0:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
92 raise Exception, 'Index file empty, there may be an error with your reference file or settings.'
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
93 except Exception, e:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
94 #clean up temp files
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
95 if os.path.exists( tmp_dir ):
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
96 shutil.rmtree( tmp_dir )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
97 stop_err( 'Error creating indexes from reference (%s), %s' % ( options.ref_file, str( e ) ) )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
98 try:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
99 # Extract all alignments from the input SAM file to BAM format ( since no region is specified, all the alignments will be extracted ).
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
100 tmp_aligns_file = tempfile.NamedTemporaryFile( dir=tmp_dir )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
101 tmp_aligns_file_name = tmp_aligns_file.name
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
102 tmp_aligns_file.close()
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
103 command = 'samtools view -bt %s -o %s %s' % ( fai_index_file_path, tmp_aligns_file_name, options.input1 )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
104 tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
105 tmp_stderr = open( tmp, 'wb' )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
106 proc = subprocess.Popen( args=command, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
107 returncode = proc.wait()
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
108 tmp_stderr.close()
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
109 # get stderr, allowing for case where it's very large
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
110 tmp_stderr = open( tmp, 'rb' )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
111 stderr = ''
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
112 buffsize = 1048576
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
113 try:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
114 while True:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
115 stderr += tmp_stderr.read( buffsize )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
116 if not stderr or len( stderr ) % buffsize != 0:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
117 break
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
118 except OverflowError:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
119 pass
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
120 tmp_stderr.close()
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
121 if returncode != 0:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
122 raise Exception, stderr
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
123 except Exception, e:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
124 #clean up temp files
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
125 if os.path.exists( tmp_dir ):
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
126 shutil.rmtree( tmp_dir )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
127 stop_err( 'Error extracting alignments from (%s), %s' % ( options.input1, str( e ) ) )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
128 try:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
129 # Sort alignments by leftmost coordinates. File <out.prefix>.bam will be created. This command
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
130 # may also create temporary files <out.prefix>.%d.bam when the whole alignment cannot be fitted
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
131 # into memory ( controlled by option -m ).
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
132 tmp_sorted_aligns_file = tempfile.NamedTemporaryFile( dir=tmp_dir )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
133 tmp_sorted_aligns_file_name = tmp_sorted_aligns_file.name
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
134 tmp_sorted_aligns_file.close()
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
135 command = 'samtools sort %s %s' % ( tmp_aligns_file_name, tmp_sorted_aligns_file_name )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
136 tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
137 tmp_stderr = open( tmp, 'wb' )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
138 proc = subprocess.Popen( args=command, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
139 returncode = proc.wait()
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
140 tmp_stderr.close()
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
141 # get stderr, allowing for case where it's very large
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
142 tmp_stderr = open( tmp, 'rb' )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
143 stderr = ''
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
144 buffsize = 1048576
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
145 try:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
146 while True:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
147 stderr += tmp_stderr.read( buffsize )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
148 if not stderr or len( stderr ) % buffsize != 0:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
149 break
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
150 except OverflowError:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
151 pass
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
152 tmp_stderr.close()
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
153 if returncode != 0:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
154 raise Exception, stderr
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
155 except Exception, e:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
156 #clean up temp files
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
157 if os.path.exists( tmp_dir ):
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
158 shutil.rmtree( tmp_dir )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
159 stop_err( 'Error sorting alignments from (%s), %s' % ( tmp_aligns_file_name, str( e ) ) )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
160 # Move tmp_aligns_file_name to our output dataset location
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
161 sorted_bam_file = '%s.bam' % tmp_sorted_aligns_file_name
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
162 shutil.move( sorted_bam_file, options.output1 )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
163 #clean up temp files
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
164 if os.path.exists( tmp_dir ):
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
165 shutil.rmtree( tmp_dir )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
166 # check that there are results in the output file
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
167 if os.path.getsize( options.output1 ) > 0:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
168 sys.stdout.write( 'SAM file converted to BAM' )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
169 else:
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
170 stop_err( 'Error creating sorted version of BAM file.' )
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
171
30fdbaccb96b Uploaded tool tarball.
devteam
parents:
diff changeset
172 if __name__=="__main__": __main__()