diff tools/samtools/bam_to_sam.py @ 0:9071e359b9a3

Uploaded
author xuebing
date Fri, 09 Mar 2012 19:37:19 -0500
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/samtools/bam_to_sam.py	Fri Mar 09 19:37:19 2012 -0500
@@ -0,0 +1,129 @@
+#!/usr/bin/env python
+"""
+Converts BAM data to sorted SAM data.
+usage: bam_to_sam.py [options]
+   --input1: SAM file to be converted
+   --output1: output dataset in bam format
+"""
+
+import optparse, os, sys, subprocess, tempfile, shutil
+from galaxy import eggs
+import pkg_resources; pkg_resources.require( "bx-python" )
+from bx.cookbook import doc_optparse
+#from galaxy import util
+
+def stop_err( msg ):
+    sys.stderr.write( '%s\n' % msg )
+    sys.exit()
+
+def __main__():
+    #Parse Command Line
+    parser = optparse.OptionParser()
+    parser.add_option( '', '--input1', dest='input1', help='The input SAM dataset' )
+    parser.add_option( '', '--output1', dest='output1', help='The output BAM dataset' )
+    parser.add_option( '', '--header', dest='header', action='store_true', default=False, help='Write SAM Header' )
+    ( options, args ) = parser.parse_args()
+
+    # output version # of tool
+    try:
+        tmp = tempfile.NamedTemporaryFile().name
+        tmp_stdout = open( tmp, 'wb' )
+        proc = subprocess.Popen( args='samtools 2>&1', shell=True, stdout=tmp_stdout )
+        tmp_stdout.close()
+        returncode = proc.wait()
+        stdout = None
+        for line in open( tmp_stdout.name, 'rb' ):
+            if line.lower().find( 'version' ) >= 0:
+                stdout = line.strip()
+                break
+        if stdout:
+            sys.stdout.write( 'Samtools %s\n' % stdout )
+        else:
+            raise Exception
+    except:
+        sys.stdout.write( 'Could not determine Samtools version\n' )
+
+    tmp_dir = tempfile.mkdtemp()
+
+    try:
+        # exit if input file empty
+        if os.path.getsize( options.input1 ) == 0:
+            raise Exception, 'Initial BAM file empty'
+        # Sort alignments by leftmost coordinates. File <out.prefix>.bam will be created. This command
+        # may also create temporary files <out.prefix>.%d.bam when the whole alignment cannot be fitted
+        # into memory ( controlled by option -m ).
+        tmp_sorted_aligns_file = tempfile.NamedTemporaryFile( dir=tmp_dir )
+        tmp_sorted_aligns_file_base = tmp_sorted_aligns_file.name
+        tmp_sorted_aligns_file_name = '%s.bam' % tmp_sorted_aligns_file.name
+        tmp_sorted_aligns_file.close()
+        command = 'samtools sort %s %s' % ( options.input1, tmp_sorted_aligns_file_base )
+        tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
+        tmp_stderr = open( tmp, 'wb' )
+        proc = subprocess.Popen( args=command, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
+        returncode = proc.wait()
+        tmp_stderr.close()
+        # get stderr, allowing for case where it's very large
+        tmp_stderr = open( tmp, 'rb' )
+        stderr = ''
+        buffsize = 1048576
+        try:
+            while True:
+                stderr += tmp_stderr.read( buffsize )
+                if not stderr or len( stderr ) % buffsize != 0:
+                    break
+        except OverflowError:
+            pass
+        tmp_stderr.close()
+        if returncode != 0:
+            raise Exception, stderr
+        # exit if sorted BAM file empty
+        if os.path.getsize( tmp_sorted_aligns_file_name) == 0:
+            raise Exception, 'Intermediate sorted BAM file empty'
+    except Exception, e:
+        #clean up temp files
+        if os.path.exists( tmp_dir ):
+            shutil.rmtree( tmp_dir )
+        stop_err( 'Error sorting alignments from (%s), %s' % ( options.input1, str( e ) ) )
+
+
+    try:
+        # Extract all alignments from the input BAM file to SAM format ( since no region is specified, all the alignments will be extracted ).
+        if options.header:
+            view_options = "-h"
+        else:
+            view_options = ""
+        command = 'samtools view %s -o %s %s' % ( view_options, options.output1, tmp_sorted_aligns_file_name )
+        tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
+        tmp_stderr = open( tmp, 'wb' )
+        proc = subprocess.Popen( args=command, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
+        returncode = proc.wait()
+        tmp_stderr.close()
+        # get stderr, allowing for case where it's very large
+        tmp_stderr = open( tmp, 'rb' )
+        stderr = ''
+        buffsize = 1048576
+        try:
+            while True:
+                stderr += tmp_stderr.read( buffsize )
+                if not stderr or len( stderr ) % buffsize != 0:
+                    break
+        except OverflowError:
+            pass
+        tmp_stderr.close()
+        if returncode != 0:
+            raise Exception, stderr
+    except Exception, e:
+        #clean up temp files
+        if os.path.exists( tmp_dir ):
+            shutil.rmtree( tmp_dir )
+        stop_err( 'Error extracting alignments from (%s), %s' % ( options.input1, str( e ) ) )
+    #clean up temp files
+    if os.path.exists( tmp_dir ):
+        shutil.rmtree( tmp_dir )
+    # check that there are results in the output file
+    if os.path.getsize( options.output1 ) > 0:
+        sys.stdout.write( 'BAM file converted to SAM' )
+    else:
+        stop_err( 'The output file is empty, there may be an error with your input file.' )
+
+if __name__=="__main__": __main__()