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