Merge pull request #3187 from mvdbeek/samtools_1_3_fix

Make sam to bam converter tool compatible with samtools >= 1.3

lib/galaxy/datatypes/converters/sam_to_bam.py

@@ -11,6 +11,7 @@ import optparse
 import tempfile
 import subprocess
 import shutil
+from distutils.version import LooseVersion
 
 CHUNK_SIZE = 2 ** 20  # 1mb
 
@@ -20,6 +21,40 @@ def cleanup_before_exit( tmp_dir ):
         shutil.rmtree( tmp_dir )
 
 
+def cmd_exists(cmd):
+    # http://stackoverflow.com/questions/5226958/which-equivalent-function-in-python
+    for path in os.environ["PATH"].split(":"):
+        if os.path.exists(os.path.join(path, cmd)):
+            return True
+    return False
+
+
+def _get_samtools_version():
+    version = '0.0.0'
+    if not cmd_exists('samtools'):
+        raise Exception('This tool needs samtools, but it is not on PATH.')
+    # Get the version of samtools via --version-only, if available
+    p = subprocess.Popen( ['samtools', '--version-only'],
+                          stdout=subprocess.PIPE,
+                          stderr=subprocess.PIPE)
+    output, error = p.communicate()
+
+    # --version-only is available
+    # Format is <version x.y.z>+htslib-<a.b.c>
+    if p.returncode == 0:
+        version = output.split('+')[0]
+        return version
+
+    output = subprocess.Popen( [ 'samtools' ], stderr=subprocess.PIPE, stdout=subprocess.PIPE ).communicate()[1]
+    lines = output.split( '\n' )
+    for line in lines:
+        if line.lower().startswith( 'version' ):
+            # Assuming line looks something like: version: 0.1.12a (r862)
+            version = line.split()[1]
+            break
+    return version
+
+
 def __main__():
     # Parse Command Line
     parser = optparse.OptionParser()
@@ -33,7 +68,7 @@ def __main__():
     # convert to SAM
     unsorted_bam_filename = os.path.join( tmp_dir, 'unsorted.bam' )
     unsorted_stderr_filename = os.path.join( tmp_dir, 'unsorted.stderr' )
-    cmd = 'samtools view -bS "%s" > "%s"' % ( input_filename, unsorted_bam_filename )
+    cmd = "samtools view -bS '%s' > '%s'" % ( input_filename, unsorted_bam_filename )
     proc = subprocess.Popen( args=cmd, stderr=open( unsorted_stderr_filename, 'wb' ), shell=True, cwd=tmp_dir )
     return_code = proc.wait()
     if return_code:
@@ -52,7 +87,12 @@ def __main__():
     # sort sam, so indexing will not fail
     sorted_stderr_filename = os.path.join( tmp_dir, 'sorted.stderr' )
     sorting_prefix = os.path.join( tmp_dir, 'sorted_bam' )
-    cmd = 'samtools sort -o "%s" "%s" > "%s"' % ( unsorted_bam_filename, sorting_prefix, output_filename )
+    # samtools changed sort command arguments (starting from version 1.3)
+    samtools_version = LooseVersion(_get_samtools_version())
+    if samtools_version < LooseVersion('1.0'):
+        cmd = "samtools sort -o '%s' '%s' > '%s'" % ( unsorted_bam_filename, sorting_prefix, output_filename )
+    else:
+        cmd = "samtools sort -T '%s' '%s' > '%s'" % ( sorting_prefix, unsorted_bam_filename, output_filename )
     proc = subprocess.Popen( args=cmd, stderr=open( sorted_stderr_filename, 'wb' ), shell=True, cwd=tmp_dir )
     return_code = proc.wait()
 

test-data/sam_with_header.sam

@@ -0,0 +1,14 @@
+@SQ	SN:ref	LN:45
+@SQ	SN:ref2	LN:40
+r001	163	ref	7	30	8M4I4M1D3M	=	37	39	TTAGATAAAGAGGATACTG	*	XX:B:S,12561,2,20,112
+r002	0	ref	9	30	1S2I6M1P1I1P1I4M2I	*	0	0	AAAAGATAAGGGATAAA	*
+r003	0	ref	9	30	5H6M	*	0	0	AGCTAA	*
+r004	0	ref	16	30	6M14N1I5M	*	0	0	ATAGCTCTCAGC	*
+r003	16	ref	29	30	6H5M	*	0	0	TAGGC	*
+r001	83	ref	37	30	9M	=	7	-39	CAGCGCCAT	*
+x1	0	ref2	1	30	20M	*	0	0	aggttttataaaacaaataa	????????????????????
+x2	0	ref2	2	30	21M	*	0	0	ggttttataaaacaaataatt	?????????????????????
+x3	0	ref2	6	30	9M4I13M	*	0	0	ttataaaacAAATaattaagtctaca	??????????????????????????
+x4	0	ref2	10	30	25M	*	0	0	CaaaTaattaagtctacagagcaac	?????????????????????????
+x5	0	ref2	12	30	24M	*	0	0	aaTaattaagtctacagagcaact	????????????????????????
+x6	0	ref2	14	30	23M	*	0	0	Taattaagtctacagagcaacta	???????????????????????

test/functional/tools/sam_to_bam.xml

@@ -0,0 +1,19 @@
+<tool id="sam_to_conversion" name="Test sam to bam conversion">
+    <command>
+        cat '$input1' > '$out_file1'
+    </command>
+    <inputs>
+        <param name="input1" type="data" format="bam" label="Concatenate Dataset"/>
+    </inputs>
+    <outputs>
+        <data name="out_file1" format="bam"/>
+    </outputs>
+    <tests>
+        <test>
+            <param name="input1" value="sam_with_header.sam" ftype="sam"/>
+            <output name="out_file1" file="bam_from_sam.bam"/>
+        </test>
+    </tests>
+    <help>
+    </help>
+</tool>