[2] | 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 |
---|
| 6 | --dbkey: dbkey value |
---|
| 7 | --ref_file: Reference file if choosing from history |
---|
| 8 | --output1: output dataset in bam format |
---|
| 9 | --index_dir: GALAXY_DATA_INDEX_DIR |
---|
| 10 | """ |
---|
| 11 | |
---|
| 12 | import optparse, os, sys, subprocess, tempfile, shutil, gzip |
---|
| 13 | from galaxy import eggs |
---|
| 14 | import pkg_resources; pkg_resources.require( "bx-python" ) |
---|
| 15 | from bx.cookbook import doc_optparse |
---|
| 16 | from galaxy import util |
---|
| 17 | |
---|
| 18 | def stop_err( msg ): |
---|
| 19 | sys.stderr.write( '%s\n' % msg ) |
---|
| 20 | sys.exit() |
---|
| 21 | |
---|
| 22 | def check_seq_file( dbkey, cached_seqs_pointer_file ): |
---|
| 23 | seq_path = '' |
---|
| 24 | for line in open( cached_seqs_pointer_file ): |
---|
| 25 | line = line.rstrip( '\r\n' ) |
---|
| 26 | if line and not line.startswith( '#' ) and line.startswith( 'index' ): |
---|
| 27 | fields = line.split( '\t' ) |
---|
| 28 | if len( fields ) < 3: |
---|
| 29 | continue |
---|
| 30 | if fields[1] == dbkey: |
---|
| 31 | seq_path = fields[2].strip() |
---|
| 32 | break |
---|
| 33 | return seq_path |
---|
| 34 | |
---|
| 35 | def __main__(): |
---|
| 36 | #Parse Command Line |
---|
| 37 | parser = optparse.OptionParser() |
---|
| 38 | parser.add_option( '', '--input1', dest='input1', help='The input SAM dataset' ) |
---|
| 39 | parser.add_option( '', '--dbkey', dest='dbkey', help='The build of the reference dataset' ) |
---|
| 40 | parser.add_option( '', '--ref_file', dest='ref_file', help='The reference dataset from the history' ) |
---|
| 41 | parser.add_option( '', '--output1', dest='output1', help='The output BAM dataset' ) |
---|
| 42 | parser.add_option( '', '--index_dir', dest='index_dir', help='GALAXY_DATA_INDEX_DIR' ) |
---|
| 43 | ( options, args ) = parser.parse_args() |
---|
| 44 | |
---|
| 45 | cached_seqs_pointer_file = '%s/sam_fa_indices.loc' % options.index_dir |
---|
| 46 | if not os.path.exists( cached_seqs_pointer_file ): |
---|
| 47 | stop_err( 'The required file (%s) does not exist.' % cached_seqs_pointer_file ) |
---|
| 48 | # If found for the dbkey, seq_path will look something like /galaxy/data/equCab2/sam_index/equCab2.fa, |
---|
| 49 | # and the equCab2.fa file will contain fasta sequences. |
---|
| 50 | seq_path = check_seq_file( options.dbkey, cached_seqs_pointer_file ) |
---|
| 51 | tmp_dir = tempfile.mkdtemp() |
---|
| 52 | if options.ref_file == 'None': |
---|
| 53 | # We're using locally cached reference sequences( e.g., /galaxy/data/equCab2/sam_index/equCab2.fa ). |
---|
| 54 | # The indexes for /galaxy/data/equCab2/sam_index/equCab2.fa will be contained in |
---|
| 55 | # a file named /galaxy/data/equCab2/sam_index/equCab2.fa.fai |
---|
| 56 | fai_index_file_base = seq_path |
---|
| 57 | fai_index_file_path = '%s.fai' % seq_path |
---|
| 58 | if not os.path.exists( fai_index_file_path ): |
---|
| 59 | #clean up temp files |
---|
| 60 | if os.path.exists( tmp_dir ): |
---|
| 61 | shutil.rmtree( tmp_dir ) |
---|
| 62 | stop_err( 'No sequences are available for build (%s), request them by reporting this error.' % options.dbkey ) |
---|
| 63 | else: |
---|
| 64 | try: |
---|
| 65 | # Create indexes for history reference ( e.g., ~/database/files/000/dataset_1.dat ) using samtools faidx, which will: |
---|
| 66 | # - index reference sequence in the FASTA format or extract subsequence from indexed reference sequence |
---|
| 67 | # - if no region is specified, faidx will index the file and create <ref.fasta>.fai on the disk |
---|
| 68 | # - if regions are specified, the subsequences will be retrieved and printed to stdout in the FASTA format |
---|
| 69 | # - the input file can be compressed in the RAZF format. |
---|
| 70 | # IMPORTANT NOTE: a real weakness here is that we are creating indexes for the history dataset |
---|
| 71 | # every time we run this tool. It would be nice if we could somehow keep track of user's specific |
---|
| 72 | # index files so they could be re-used. |
---|
| 73 | fai_index_file_base = tempfile.NamedTemporaryFile( dir=tmp_dir ).name |
---|
| 74 | # At this point, fai_index_file_path will look something like /tmp/dataset_13.dat |
---|
| 75 | os.symlink( options.ref_file, fai_index_file_base ) |
---|
| 76 | fai_index_file_path = '%s.fai' % fai_index_file_base |
---|
| 77 | command = 'samtools faidx %s' % fai_index_file_base |
---|
| 78 | tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name |
---|
| 79 | tmp_stderr = open( tmp, 'wb' ) |
---|
| 80 | proc = subprocess.Popen( args=command, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() ) |
---|
| 81 | returncode = proc.wait() |
---|
| 82 | tmp_stderr.close() |
---|
| 83 | # get stderr, allowing for case where it's very large |
---|
| 84 | tmp_stderr = open( tmp, 'rb' ) |
---|
| 85 | stderr = '' |
---|
| 86 | buffsize = 1048576 |
---|
| 87 | try: |
---|
| 88 | while True: |
---|
| 89 | stderr += tmp_stderr.read( buffsize ) |
---|
| 90 | if not stderr or len( stderr ) % buffsize != 0: |
---|
| 91 | break |
---|
| 92 | except OverflowError: |
---|
| 93 | pass |
---|
| 94 | tmp_stderr.close() |
---|
| 95 | if returncode != 0: |
---|
| 96 | raise Exception, stderr |
---|
| 97 | if os.path.getsize( fai_index_file_path ) == 0: |
---|
| 98 | raise Exception, 'Index file empty, there may be an error with your reference file or settings.' |
---|
| 99 | except Exception, e: |
---|
| 100 | #clean up temp files |
---|
| 101 | if os.path.exists( tmp_dir ): |
---|
| 102 | shutil.rmtree( tmp_dir ) |
---|
| 103 | stop_err( 'Error creating indexes from reference (%s), %s' % ( options.ref_file, str( e ) ) ) |
---|
| 104 | try: |
---|
| 105 | # Extract all alignments from the input SAM file to BAM format ( since no region is specified, all the alignments will be extracted ). |
---|
| 106 | tmp_aligns_file = tempfile.NamedTemporaryFile( dir=tmp_dir ) |
---|
| 107 | tmp_aligns_file_name = tmp_aligns_file.name |
---|
| 108 | tmp_aligns_file.close() |
---|
| 109 | # IMPORTANT NOTE: for some reason the samtools view command gzips the resulting bam file without warning, |
---|
| 110 | # and the docs do not currently state that this occurs ( very bad ). |
---|
| 111 | command = 'samtools view -bt %s -o %s %s' % ( fai_index_file_path, tmp_aligns_file_name, options.input1 ) |
---|
| 112 | tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name |
---|
| 113 | tmp_stderr = open( tmp, 'wb' ) |
---|
| 114 | proc = subprocess.Popen( args=command, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() ) |
---|
| 115 | returncode = proc.wait() |
---|
| 116 | tmp_stderr.close() |
---|
| 117 | # get stderr, allowing for case where it's very large |
---|
| 118 | tmp_stderr = open( tmp, 'rb' ) |
---|
| 119 | stderr = '' |
---|
| 120 | buffsize = 1048576 |
---|
| 121 | try: |
---|
| 122 | while True: |
---|
| 123 | stderr += tmp_stderr.read( buffsize ) |
---|
| 124 | if not stderr or len( stderr ) % buffsize != 0: |
---|
| 125 | break |
---|
| 126 | except OverflowError: |
---|
| 127 | pass |
---|
| 128 | tmp_stderr.close() |
---|
| 129 | if returncode != 0: |
---|
| 130 | raise Exception, stderr |
---|
| 131 | if len( open( tmp_aligns_file_name ).read() ) == 0: |
---|
| 132 | raise Exception, 'Initial BAM file empty' |
---|
| 133 | except Exception, e: |
---|
| 134 | #clean up temp files |
---|
| 135 | if os.path.exists( tmp_dir ): |
---|
| 136 | shutil.rmtree( tmp_dir ) |
---|
| 137 | stop_err( 'Error extracting alignments from (%s), %s' % ( options.input1, str( e ) ) ) |
---|
| 138 | try: |
---|
| 139 | # Sort alignments by leftmost coordinates. File <out.prefix>.bam will be created. This command |
---|
| 140 | # may also create temporary files <out.prefix>.%d.bam when the whole alignment cannot be fitted |
---|
| 141 | # into memory ( controlled by option -m ). |
---|
| 142 | tmp_sorted_aligns_file = tempfile.NamedTemporaryFile( dir=tmp_dir ) |
---|
| 143 | tmp_sorted_aligns_file_name = tmp_sorted_aligns_file.name |
---|
| 144 | tmp_sorted_aligns_file.close() |
---|
| 145 | command = 'samtools sort %s %s' % ( tmp_aligns_file_name, tmp_sorted_aligns_file_name ) |
---|
| 146 | tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name |
---|
| 147 | tmp_stderr = open( tmp, 'wb' ) |
---|
| 148 | proc = subprocess.Popen( args=command, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() ) |
---|
| 149 | returncode = proc.wait() |
---|
| 150 | tmp_stderr.close() |
---|
| 151 | # get stderr, allowing for case where it's very large |
---|
| 152 | tmp_stderr = open( tmp, 'rb' ) |
---|
| 153 | stderr = '' |
---|
| 154 | buffsize = 1048576 |
---|
| 155 | try: |
---|
| 156 | while True: |
---|
| 157 | stderr += tmp_stderr.read( buffsize ) |
---|
| 158 | if not stderr or len( stderr ) % buffsize != 0: |
---|
| 159 | break |
---|
| 160 | except OverflowError: |
---|
| 161 | pass |
---|
| 162 | tmp_stderr.close() |
---|
| 163 | if returncode != 0: |
---|
| 164 | raise Exception, stderr |
---|
| 165 | except Exception, e: |
---|
| 166 | #clean up temp files |
---|
| 167 | if os.path.exists( tmp_dir ): |
---|
| 168 | shutil.rmtree( tmp_dir ) |
---|
| 169 | stop_err( 'Error sorting alignments from (%s), %s' % ( tmp_aligns_file_name, str( e ) ) ) |
---|
| 170 | # Move tmp_aligns_file_name to our output dataset location |
---|
| 171 | sorted_bam_file = '%s.bam' % tmp_sorted_aligns_file_name |
---|
| 172 | if os.path.getsize( sorted_bam_file ) == 0: |
---|
| 173 | #clean up temp files |
---|
| 174 | if os.path.exists( tmp_dir ): |
---|
| 175 | shutil.rmtree( tmp_dir ) |
---|
| 176 | stop_err( 'Error creating sorted version of BAM file' ) |
---|
| 177 | shutil.move( sorted_bam_file, options.output1 ) |
---|
| 178 | #clean up temp files |
---|
| 179 | if os.path.exists( tmp_dir ): |
---|
| 180 | shutil.rmtree( tmp_dir ) |
---|
| 181 | # check that there are results in the output file |
---|
| 182 | if os.path.getsize( options.output1 ) > 0: |
---|
| 183 | sys.stdout.write( 'SAM file converted to BAM' ) |
---|
| 184 | else: |
---|
| 185 | stop_err( 'The output file is empty, there may be an error with your input file or settings.' ) |
---|
| 186 | |
---|
| 187 | if __name__=="__main__": __main__() |
---|