/tools/samtools/sam_to_bam.py

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