diff --git a/tools/samtools/sam_merge.py b/tools/samtools/sam_merge.py
index 6f6c7db926c..d62e9e7f97d 100644
--- a/tools/samtools/sam_merge.py
+++ b/tools/samtools/sam_merge.py
@@ -1,21 +1,38 @@
#! /usr/bin/python
-import os, sys
+"""
+Merges any number of BAM files
+usage: %prog [options]
+ input1
+ output1
+ input2
+ [input3[,input4[,input5[,...]]]]
+"""
+
+import os, subprocess, sys
def stop_err( msg ):
- sys.stderr.write( msg )
+ sys.stderr.write( '%s\n' % msg )
sys.exit()
def __main__():
infile = sys.argv[1]
outfile = sys.argv[2]
if len( sys.argv ) < 3:
- stop_err( 'No files to merge' )
+ stop_err( 'There are not enough files to merge' )
filenames = sys.argv[3:]
- cmd1 = 'samtools merge %s %s %s' % (outfile, infile, ' '.join(filenames))
+ cmd = 'samtools merge %s %s %s' % ( outfile, infile, ' '.join( filenames ) )
try:
- os.system(cmd1)
- except Exception, eq:
- stop_err('Error running SAMtools merge tool\n' + str(eq))
-
-if __name__ == "__main__" : __main__()
\ No newline at end of file
+ proc = subprocess.Popen( args=cmd, shell=True, stderr=subprocess.PIPE )
+ returncode = proc.wait()
+ stderr = proc.stderr.read()
+ if returncode != 0:
+ raise Exception, stderr
+ except Exception, e:
+ stop_err( 'Error running SAMtools merge tool\n' + str( e ) )
+ if os.path.getsize( outfile ) > 0:
+ sys.stdout.write( '%s files merged.' % ( len( sys.argv ) - 2 ) )
+ else:
+ stop_err( 'The output file is empty, there may be an error with one of your input files.' )
+
+if __name__ == "__main__" : __main__()
diff --git a/tools/samtools/sam_merge.xml b/tools/samtools/sam_merge.xml
index 621ef0e6641..b18dd681291 100644
--- a/tools/samtools/sam_merge.xml
+++ b/tools/samtools/sam_merge.xml
@@ -19,7 +19,48 @@
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
**What it does**
diff --git a/tools/samtools/sam_pileup.py b/tools/samtools/sam_pileup.py
index a03487ff2a4..07baa1b83d1 100644
--- a/tools/samtools/sam_pileup.py
+++ b/tools/samtools/sam_pileup.py
@@ -22,80 +22,94 @@ usage: %prog [options]
"""
-import os, sys, tempfile
+import os, shutil, subprocess, sys, tempfile
from galaxy import eggs
import pkg_resources; pkg_resources.require( "bx-python" )
from bx.cookbook import doc_optparse
def stop_err( msg ):
- sys.stderr.write( msg )
+ sys.stderr.write( '%s\n' % msg )
sys.exit()
def check_seq_file( dbkey, GALAXY_DATA_INDEX_DIR ):
- seq_file = "%s/sam_fa_indices.loc" % GALAXY_DATA_INDEX_DIR
- seq_path = ''
- for line in open( seq_file ):
+ seqFile = '%s/sam_fa_indices.loc' % GALAXY_DATA_INDEX_DIR
+ seqPath = ''
+ for line in open( seqFile ):
line = line.rstrip( '\r\n' )
- if line and not line.startswith( "#" ) and line.startswith( 'index' ):
+ if line and not line.startswith( '#' ) and line.startswith( 'index' ):
fields = line.split( '\t' )
if len( fields ) < 3:
continue
if fields[1] == dbkey:
- seq_path = fields[2].strip()
+ seqPath = fields[2].strip()
break
- return seq_path
-
+ return seqPath
+
def __main__():
#Parse Command Line
options, args = doc_optparse.parse( __doc__ )
- seq_path = check_seq_file( options.dbkey, options.indexDir )
- tmp_dir = tempfile.gettempdir()
- os.chdir(tmp_dir)
- tmpf0 = tempfile.NamedTemporaryFile(dir=tmp_dir)
- tmpf0bam = '%s.bam' % tmpf0.name
- tmpf0bambai = '%s.bam.bai' % tmpf0.name
- tmpf1 = tempfile.NamedTemporaryFile(dir=tmp_dir)
- tmpf1fai = '%s.fai' % tmpf1.name
- opts = '%s %s -M %s' % (('','-s')[options.lastCol=='yes'], ('','-i')[options.indels=='yes'], options.mapCap)
- if options.consensus == 'yes':
- opts += ' -c -T %s -N %s -r %s -I %s' % (options.theta, options.hapNum, options.fraction, options.phredProb)
- cmd1 = None
- cmd2 = 'cp %s %s; cp %s %s' % (options.input1, tmpf0bam, options.bamIndex, tmpf0bambai)
- cmd3 = 'samtools pileup %s -f %s %s > %s 2> /dev/null'
- if options.ref =='indexed':
- full_path = "%s.fai" % seq_path
- if not os.path.exists( full_path ):
- stop_err( "No sequences are available for '%s', request them by reporting this error." % options.dbkey )
- cmd3 = cmd3 % (opts, seq_path, tmpf0bam, options.output1)
- elif options.ref == 'history':
- cmd1 = 'cp %s %s; samtools faidx %s' % (options.ownFile, tmpf1.name, tmpf1.name)
- cmd3 = cmd3 % (opts, tmpf1.name, tmpf0bam, options.output1)
- # index reference if necessary
- if cmd1:
- try:
- os.system(cmd1)
- if options.ref == 'history' and not os.path.exists( tmpf1fai ):
- stop_err( "Problem creating index file from history item." )
- except Exception, eq:
- stop_err('Error handling reference sequence\n' + str(eq))
- # copy bam index to working directory
- try:
- os.system(cmd2)
- except Exception, eq:
- stop_err('Error moving files to temp directory\n' + str(eq))
- # perform pileup command
- try:
- os.system(cmd3)
- except Exception, eq:
- stop_err('Error running SAMtools pileup tool\n' + str(eq))
- # clean up temp files
- tmpf1.close()
+ seqPath = check_seq_file( options.dbkey, options.indexDir )
+ #prepare file names
+ tmpDir = tempfile.mkdtemp()
+ tmpf0 = tempfile.NamedTemporaryFile( dir=tmpDir )
+ tmpf0_name = tmpf0.name
tmpf0.close()
- if os.path.exists(tmpf0bam):
- os.remove(tmpf0bam)
- if os.path.exists(tmpf0bambai):
- os.remove(tmpf0bambai)
- if os.path.exists(tmpf1fai):
- os.remove(tmpf1fai)
-
+ tmpf0bam_name = '%s.bam' % tmpf0_name
+ tmpf0bambai_name = '%s.bam.bai' % tmpf0_name
+ tmpf1 = tempfile.NamedTemporaryFile( dir=tmpDir )
+ tmpf1_name = tmpf1.name
+ tmpf1.close()
+ tmpf1fai_name = '%s.fai' % tmpf1_name
+ #link bam and bam index to working directory (can't move because need to leave original)
+ os.symlink( options.input1, tmpf0bam_name )
+ os.symlink( options.bamIndex, tmpf0bambai_name )
+ #get parameters for pileup command
+ if options.lastCol == 'yes':
+ lastCol = '-s'
+ else:
+ lastCol = ''
+ if options.indels == 'yes':
+ indels = '-i'
+ else:
+ indels = ''
+ opts = '%s %s -M %s' % ( lastCol, indels, options.mapCap )
+ if options.consensus == 'yes':
+ opts += ' -c -T %s -N %s -r %s -I %s' % ( options.theta, options.hapNum, options.fraction, options.phredProb )
+ #prepare basic pileup command
+ cmd = 'samtools pileup %s -f %s %s > %s'
+ try:
+ #index reference if necessary and prepare pileup command
+ if options.ref == 'indexed':
+ if not os.path.exists( "%s.fai" % seqPath ):
+ raise Exception, "No sequences are available for '%s', request them by reporting this error." % options.dbkey
+ cmd = cmd % ( opts, seqPath, tmpf0bam_name, options.output1 )
+ elif options.ref == 'history':
+ os.symlink( options.ownFile, tmpf1_name )
+ cmdIndex = 'samtools faidx %s' % ( tmpf1_name )
+ proc = subprocess.Popen( args=cmdIndex, shell=True, cwd=tmpDir, stderr=subprocess.PIPE )
+ returncode = proc.wait()
+ stderr = proc.stderr.read()
+ #did index succeed?
+ if returncode != 0:
+ raise Exception, 'Error creating index file\n' + stderr
+ cmd = cmd % ( opts, tmpf1_name, tmpf0bam_name, options.output1 )
+ #perform pileup command
+ proc = subprocess.Popen( args=cmd, shell=True, cwd=tmpDir, stderr=subprocess.PIPE )
+ returncode = proc.wait()
+ #did it succeed?
+ stderr = proc.stderr.read()
+ if returncode != 0:
+ raise Exception, stderr
+ except Exception, e:
+ stop_err( 'Error running Samtools pileup tool\n' + str( e ) )
+ finally:
+ #clean up temp files
+ if os.path.exists( tmpDir ):
+ shutil.rmtree( tmpDir )
+ # check that there are results in the output file
+ if os.path.getsize( options.output1 ) > 0:
+ sys.stdout.write( 'Converted BAM to pileup' )
+ else:
+ stop_err( 'The output file is empty. Your input file may have had no matches, or there may be an error with your input file or settings.' )
+
if __name__ == "__main__" : __main__()
diff --git a/tools/samtools/sam_pileup.xml b/tools/samtools/sam_pileup.xml
index 3fce5731bb3..f2a8362f215 100644
--- a/tools/samtools/sam_pileup.xml
+++ b/tools/samtools/sam_pileup.xml
@@ -1,16 +1,16 @@
from BAM dataset
- sam_pileup.py
- --input1=$input1
- --output=$output1
- --ref=$refOrHistory.reference
- #if $refOrHistory.reference == "history":
- --ownFile=$refOrHistory.ownFile
- #else:
- --ownFile="None"
- #end if
- --dbkey=${input1.metadata.dbkey}
+ sam_pileup.py
+ --input1=$input1
+ --output=$output1
+ --ref=$refOrHistory.reference
+ #if $refOrHistory.reference == "history":
+ --ownFile=$refOrHistory.ownFile
+ #else:
+ --ownFile="None"
+ #end if
+ --dbkey=${input1.metadata.dbkey}
--indexDir=${GALAXY_DATA_INDEX_DIR}
--bamIndex=${input1.metadata.bam_index}
--lastCol=$lastCol
@@ -40,12 +40,12 @@
-
+
-
+
@@ -72,10 +72,45 @@
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-**What it does**
+
+**What it does**
Uses SAMTools_' pileup command to produce a pileup dataset from a provided BAM dataset. It generates two types of pileup datasets depending on the specified options. If *Call consensus according to MAQ model?* option is set to **No**, the tool produces simple pileup. If the option is set to **Yes**, a ten column pileup dataset with consensus is generated. Both types of datasets are briefly summarized below.
@@ -92,7 +127,7 @@ The description of pileup format below is largely based on information that can
**Six column pileup**::
1 2 3 4 5 6
- ---------------------------------
+ ---------------------------------
chrM 412 A 2 ., II
chrM 413 G 4 ..t, IIIH
chrM 414 C 4 ...a III2
@@ -124,7 +159,7 @@ The `ten-column` (consensus_) pileup incorporates additional consensus informati
where::
Column Definition
- ------- ----------------------------
+ ------- --------------------------------------------------------
1 Chromosome
2 Position (1-based)
3 Reference base at that position
diff --git a/tools/samtools/sam_to_bam.py b/tools/samtools/sam_to_bam.py
index 114d1b76b07..f11dba69e0c 100644
--- a/tools/samtools/sam_to_bam.py
+++ b/tools/samtools/sam_to_bam.py
@@ -16,14 +16,14 @@ from bx.cookbook import doc_optparse
from galaxy import util
def stop_err( msg ):
- sys.stderr.write( "%s\n" % msg )
+ sys.stderr.write( '%s\n' % msg )
sys.exit()
def check_seq_file( dbkey, cached_seqs_pointer_file ):
seq_path = ''
for line in open( cached_seqs_pointer_file ):
line = line.rstrip( '\r\n' )
- if line and not line.startswith( "#" ) and line.startswith( 'index' ):
+ if line and not line.startswith( '#' ) and line.startswith( 'index' ):
fields = line.split( '\t' )
if len( fields ) < 3:
continue
@@ -31,7 +31,7 @@ def check_seq_file( dbkey, cached_seqs_pointer_file ):
seq_path = fields[2].strip()
break
return seq_path
-
+
def __main__():
#Parse Command Line
parser = optparse.OptionParser()
@@ -42,20 +42,24 @@ def __main__():
parser.add_option( '', '--index_dir', dest='index_dir', help='GALAXY_DATA_INDEX_DIR' )
( options, args ) = parser.parse_args()
- cached_seqs_pointer_file = "%s/sam_fa_indices.loc" % options.index_dir
+ cached_seqs_pointer_file = '%s/sam_fa_indices.loc' % options.index_dir
if not os.path.exists( cached_seqs_pointer_file ):
- stop_err( "The required file (%s) does not exist." % cached_seqs_pointer_file )
+ stop_err( 'The required file (%s) does not exist.' % cached_seqs_pointer_file )
# If found for the dbkey, seq_path will look something like /depot/data2/galaxy/equCab2/sam_index/equCab2.fa,
# and the equCab2.fa file will contain fasta sequences.
seq_path = check_seq_file( options.dbkey, cached_seqs_pointer_file )
- tmp_dir = tempfile.gettempdir()
- if options.ref_file == "None":
+ tmp_dir = tempfile.mkdtemp()
+ if options.ref_file == 'None':
# We're using locally cached reference sequences( e.g., /depot/data2/galaxy/equCab2/sam_index/equCab2.fa ).
# The indexes for /depot/data2/galaxy/equCab2/sam_index/equCab2.fa will be contained in
# a file named /depot/data2/galaxy/equCab2/sam_index/equCab2.fa.fai
- fai_index_file_path = "%s.fai" % seq_path
+ fai_index_file_base = seq_path
+ fai_index_file_path = '%s.fai' % seq_path
if not os.path.exists( fai_index_file_path ):
- stop_err( "No sequences are available for build (%s), request them by reporting this error." % options.dbkey )
+ #clean up temp files
+ if os.path.exists( tmp_dir ):
+ shutil.rmtree( tmp_dir )
+ stop_err( 'No sequences are available for build (%s), request them by reporting this error.' % options.dbkey )
else:
try:
# Create indexes for history reference ( e.g., ~/database/files/000/dataset_1.dat ) using samtools faidx, which will:
@@ -66,48 +70,76 @@ def __main__():
# IMPORTANT NOTE: a real weakness here is that we are creating indexes for the history dataset
# every time we run this tool. It would be nice if we could somehow keep track of user's specific
# index files so they could be re-used.
- fai_index_file_path = os.path.join( tmp_dir, os.path.basename( options.ref_file ) )
+ fai_index_file_base = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
# At this point, fai_index_file_path will look something like /tmp/dataset_13.dat
- os.symlink( options.ref_file, fai_index_file_path )
- command = "samtools faidx %s 2>/dev/null" % fai_index_file_path
- proc = subprocess.Popen( args=command, shell=True )
- proc.wait()
+ os.symlink( options.ref_file, fai_index_file_base )
+ fai_index_file_path = '%s.fai' % fai_index_file_base
+ command = 'samtools faidx %s' % fai_index_file_base
+ proc = subprocess.Popen( args=command, shell=True, cwd=tmp_dir, stderr=subprocess.PIPE )
+ returncode = proc.wait()
+ stderr = proc.stderr.read()
+ if returncode != 0:
+ raise Exception, stderr
+ if len( open( fai_index_file_path ).read().strip() ) == 0:
+ raise Exception, 'Index file empty, there may be an error with your reference file or settings.'
except Exception, e:
+ #clean up temp files
+ if os.path.exists( tmp_dir ):
+ shutil.rmtree( tmp_dir )
stop_err( 'Error creating indexes from reference (%s), %s' % ( options.ref_file, str( e ) ) )
try:
# Extract all alignments from the input SAM file to BAM format ( since no region is specified, all the alignments will be extracted ).
- tmp_aligns_file = tempfile.NamedTemporaryFile()
+ tmp_aligns_file = tempfile.NamedTemporaryFile( dir=tmp_dir )
tmp_aligns_file_name = tmp_aligns_file.name
tmp_aligns_file.close()
# IMPORTANT NOTE: for some reason the samtools view command gzips the resulting bam file without warning,
# and the docs do not currently state that this occurs ( very bad ).
- command = "samtools view -bt %s -o %s %s 2>/dev/null" % ( fai_index_file_path, tmp_aligns_file_name, options.input1 )
- proc = subprocess.Popen( args=command, shell=True )
- proc.wait()
+ command = 'samtools view -bt %s -o %s %s' % ( fai_index_file_path, tmp_aligns_file_name, options.input1 )
+ proc = subprocess.Popen( args=command, shell=True, cwd=tmp_dir, stderr=subprocess.PIPE )
+ returncode = proc.wait()
+ stderr = proc.stderr.read()
+ if returncode != 0:
+ raise Exception, stderr
+ if len( open( tmp_aligns_file_name ).read() ) == 0:
+ raise Exception, 'Initial BAM file empty'
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 ) ) )
try:
# Sort alignments by leftmost coordinates. File .bam will be created. This command
# may also create temporary files .%d.bam when the whole alignment cannot be fitted
# into memory ( controlled by option -m ).
- tmp_sorted_aligns_file = tempfile.NamedTemporaryFile()
+ tmp_sorted_aligns_file = tempfile.NamedTemporaryFile( dir=tmp_dir )
tmp_sorted_aligns_file_name = tmp_sorted_aligns_file.name
tmp_sorted_aligns_file.close()
- command = "samtools sort %s %s 2>/dev/null" % ( tmp_aligns_file_name, tmp_sorted_aligns_file_name )
- proc = subprocess.Popen( args=command, shell=True )
- proc.wait()
+ command = 'samtools sort %s %s' % ( tmp_aligns_file_name, tmp_sorted_aligns_file_name )
+ proc = subprocess.Popen( args=command, shell=True, cwd=tmp_dir, stderr=subprocess.PIPE )
+ returncode = proc.wait()
+ stderr = proc.stderr.read()
+ 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 sorting alignments from (%s), %s' % ( tmp_aligns_file_name, str( e ) ) )
# Move tmp_aligns_file_name to our output dataset location
sorted_bam_file = '%s.bam' % tmp_sorted_aligns_file_name
+ if os.path.getsize( sorted_bam_file ) == 0:
+ #clean up temp files
+ if os.path.exists( tmp_dir ):
+ shutil.rmtree( tmp_dir )
+ stop_err( 'Error creating sorted version of BAM file' )
shutil.move( sorted_bam_file, options.output1 )
- if options.ref_file != "None":
- # Remove the symlink from /tmp/dataset_13.dat to ~/database/files/000/dataset_13.dat
- os.unlink( fai_index_file_path )
- # Remove the index file
- index_file_name = '%s.fai' % fai_index_file_path
- os.unlink( index_file_name )
- # Remove the tmp_aligns_file_name
- os.unlink( tmp_aligns_file_name )
+ #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( 'SAM file converted to BAM' )
+ else:
+ stop_err( 'The output file is empty, there may be an error with your input file or settings.' )
if __name__=="__main__": __main__()
diff --git a/tools/samtools/sam_to_bam.xml b/tools/samtools/sam_to_bam.xml
index b6eea3b1675..024b880e2dd 100644
--- a/tools/samtools/sam_to_bam.xml
+++ b/tools/samtools/sam_to_bam.xml
@@ -32,11 +32,30 @@ sam_to_bam.py --input1=$source.input1 --dbkey=${input1.metadata.dbkey}
+
+
+
+
+
+
+
diff --git a/tools/sr_mapping/bwa_wrapper.py b/tools/sr_mapping/bwa_wrapper.py
index 57198586327..f27eba3a682 100644
--- a/tools/sr_mapping/bwa_wrapper.py
+++ b/tools/sr_mapping/bwa_wrapper.py
@@ -3,7 +3,7 @@
"""
Runs BWA on single-end or paired-end data.
Produces a SAM file containing the mappings.
-Works with BWA version 0.5.3.
+Works with BWA version 0.5.3-0.5.5.
usage: bwa_wrapper.py [options]
-t, --threads=t: The number of threads to use
@@ -34,12 +34,12 @@ usage: bwa_wrapper.py [options]
-H, --suppressHeader=h: Suppress header
"""
-import optparse, os, shutil, sys, tempfile
+import optparse, os, shutil, subprocess, sys, tempfile
def stop_err( msg ):
- sys.stderr.write( "%s\n" % msg )
+ sys.stderr.write( '%s\n' % msg )
sys.exit()
-
+
def __main__():
#Parse Command Line
parser = optparse.OptionParser()
@@ -70,14 +70,16 @@ def __main__():
parser.add_option( '-D', '--dbkey', dest='dbkey', help='Dbkey for reference genome' )
parser.add_option( '-H', '--suppressHeader', dest='suppressHeader', help='Suppress header' )
(options, args) = parser.parse_args()
- # make temp directory for placement of indices and copy reference file there
+ # make temp directory for placement of indices
tmp_index_dir = tempfile.mkdtemp()
+ tmp_dir = tempfile.mkdtemp()
# index if necessary
if options.fileSource == 'history':
- try:
- shutil.copy( options.ref, tmp_index_dir )
- except Exception, e:
- stop_err( 'Error creating temp directory for indexing purposes\n' + str( e ) )
+ ref_file = tempfile.NamedTemporaryFile( dir=tmp_index_dir )
+ ref_file_name = ref_file.name
+ ref_file.close()
+ os.symlink( options.ref, ref_file_name )
+ # determine which indexing algorithm to use, based on size
try:
size = os.stat( options.ref ).st_size
if size <= 2**30:
@@ -87,13 +89,22 @@ def __main__():
except:
indexingAlg = 'is'
indexing_cmds = '-a %s' % indexingAlg
- options.ref = os.path.join( tmp_index_dir, os.path.split( options.ref )[1] )
- cmd1 = 'bwa index %s %s 2> /dev/null' % ( indexing_cmds, options.ref )
+ cmd1 = 'bwa index %s %s' % ( indexing_cmds, ref_file_name )
try:
- os.chdir( tmp_index_dir )
- os.system( cmd1 )
+ proc = subprocess.Popen( args=cmd1, shell=True, cwd=tmp_index_dir, stderr=subprocess.PIPE )
+ returncode = proc.wait()
+ stderr = proc.stderr.read()
+ if returncode != 0:
+ raise Exception, stderr
except Exception, e:
- stop_err( 'Error indexing reference sequence\n' + str( e ) )
+ # clean up temp dirs
+ if os.path.exists( tmp_index_dir ):
+ shutil.rmtree( tmp_index_dir )
+ if os.path.exists( tmp_dir ):
+ shutil.rmtree( tmp_dir )
+ stop_err( 'Error indexing reference sequence. ' + str( e ) )
+ else:
+ ref_file_name = options.ref
# set up aligning and generate aligning command options
if options.params == 'pre_set':
aligning_cmds = '-t %s' % options.threads
@@ -116,60 +127,85 @@ def __main__():
else:
noIterSearch = ''
aligning_cmds = '-n %s -o %s -e %s -d %s -i %s %s -k %s -t %s -M %s -O %s -E %s %s %s' % \
- ( editDist, options.maxGapOpens, options.maxGapExtens, options.disallowLongDel,
- options.disallowIndel, seed, options.maxEditDistSeed, options.threads,
- options.mismatchPenalty, options.gapOpenPenalty, options.gapExtensPenalty,
+ ( editDist, options.maxGapOpens, options.maxGapExtens, options.disallowLongDel,
+ options.disallowIndel, seed, options.maxEditDistSeed, options.threads,
+ options.mismatchPenalty, options.gapOpenPenalty, options.gapExtensPenalty,
suboptAlign, noIterSearch )
if options.genAlignType == 'single':
gen_alignment_cmds = '-n %s' % options.outputTopN
elif options.genAlignType == 'paired':
gen_alignment_cmds = '-a %s -o %s' % ( options.maxInsertSize, options.maxOccurPairing )
-# print 'options.genAlignType: %s and commands: %s' % (options.genAlignType, gen_alignment_cmds)
# set up output files
- tmp_align_out = tempfile.NamedTemporaryFile()
- tmp_align_out2 = tempfile.NamedTemporaryFile()
+ tmp_align_out = tempfile.NamedTemporaryFile( dir=tmp_dir )
+ tmp_align_out_name = tmp_align_out.name
+ tmp_align_out.close()
+ tmp_align_out2 = tempfile.NamedTemporaryFile( dir=tmp_dir )
+ tmp_align_out2_name = tmp_align_out2.name
+ tmp_align_out2.close()
# prepare actual aligning and generate aligning commands
- cmd2 = 'bwa aln %s %s %s > %s 2> /dev/null' % ( aligning_cmds, options.ref, options.fastq, tmp_align_out.name )
+ cmd2 = 'bwa aln %s %s %s > %s' % ( aligning_cmds, ref_file_name, options.fastq, tmp_align_out_name )
cmd2b = ''
if options.genAlignType == 'paired':
- cmd2b = 'bwa aln %s %s %s > %s 2> /dev/null' % ( aligning_cmds, options.ref, options.rfastq, tmp_align_out2.name )
- cmd3 = 'bwa sampe %s %s %s %s %s %s >> %s 2> /dev/null' % ( gen_alignment_cmds, options.ref, tmp_align_out.name, tmp_align_out2.name, options.fastq, options.rfastq, options.output )
+ cmd2b = 'bwa aln %s %s %s > %s' % ( aligning_cmds, ref_file_name, options.rfastq, tmp_align_out2_name )
+ cmd3 = 'bwa sampe %s %s %s %s %s %s >> %s' % ( gen_alignment_cmds, ref_file_name, tmp_align_out_name, tmp_align_out2_name, options.fastq, options.rfastq, options.output )
else:
- cmd3 = 'bwa samse %s %s %s %s >> %s 2> /dev/null' % ( gen_alignment_cmds, options.ref, tmp_align_out.name, options.fastq, options.output )
- # align
+ cmd3 = 'bwa samse %s %s %s %s >> %s' % ( gen_alignment_cmds, ref_file_name, tmp_align_out_name, options.fastq, options.output )
+ # perform alignments
try:
- os.system( cmd2 )
- except Exception, e:
- stop_err( 'Error aligning sequence\n' + str( e ) )
- # and again if paired data
- try:
- if cmd2b:
- os.system( cmd2b )
- except Exception, erf:
- stop_err( 'Error aligning second sequence\n' + str( e ) )
- # generate align
- try:
- os.system( cmd3 )
- except Exception, e:
- stop_err( 'Error sequence aligning sequence\n' + str( e ) )
- # clean up temp files
- tmp_align_out.close()
- tmp_align_out2.close()
- # remove header if necessary
- if options.suppressHeader == 'true':
- tmp_out = tempfile.NamedTemporaryFile()
+ # align
try:
- shutil.move( options.output, tmp_out.name )
+ proc = subprocess.Popen( args=cmd2, shell=True, cwd=tmp_dir, stderr=subprocess.PIPE )
+ returncode = proc.wait()
+ stderr = proc.stderr.read()
+ if returncode != 0:
+ raise Exception, stderr
except Exception, e:
- stop_err( 'Error moving output file before removing headers\n' + str( e ) )
- fout = file( options.output, 'w' )
- for line in file( tmp_out.name, 'r' ):
- if not ( line.startswith( '@HD' ) or line.startswith( '@SQ' ) or line.startswith( '@RG' ) or line.startswith( '@PG' ) or line.startswith( '@CO' ) ):
- fout.write( line )
- fout.close()
- tmp_out.close()
- # clean up temp dir
- if os.path.exists( tmp_index_dir ):
- shutil.rmtree( tmp_index_dir )
-
+ raise Exception, 'Error aligning sequence. ' + str( e )
+ # and again if paired data
+ try:
+ if cmd2b:
+ proc = subprocess.Popen( args=cmd2b, shell=True, cwd=tmp_dir, stderr=subprocess.PIPE )
+ returncode = proc.wait()
+ stderr = proc.stderr.read()
+ if returncode != 0:
+ raise Exception, stderr
+ except Exception, e:
+ raise Exception, 'Error aligning second sequence. ' + str( e )
+ # generate align
+ try:
+ proc = subprocess.Popen( args=cmd3, shell=True, cwd=tmp_dir, stderr=subprocess.PIPE )
+ returncode = proc.wait()
+ stderr = proc.stderr.read()
+ if returncode != 0:
+ raise Exception, stderr
+ except Exception, e:
+ raise Exception, 'Error generating alignments. ' + str( e )
+ # remove header if necessary
+ if options.suppressHeader == 'true':
+ tmp_out = tempfile.NamedTemporaryFile( dir=tmp_dir)
+ tmp_out_name = tmp_out.name
+ tmp_out.close()
+ try:
+ shutil.move( options.output, tmp_out_name )
+ except Exception, e:
+ raise Exception, 'Error moving output file before removing headers. ' + str( e )
+ fout = file( options.output, 'w' )
+ for line in file( tmp_out.name, 'r' ):
+ if not ( line.startswith( '@HD' ) or line.startswith( '@SQ' ) or line.startswith( '@RG' ) or line.startswith( '@PG' ) or line.startswith( '@CO' ) ):
+ fout.write( line )
+ fout.close()
+ # check that there are results in the output file
+ if os.path.getsize( options.output ) > 0:
+ sys.stdout.write( 'BWA run on %s-end data' % options.genAlignType )
+ else:
+ raise Exception, 'The output file is empty. You may simply have no matches, or there may be an error with your input file or settings.'
+ except Exception, e:
+ stop_err( 'The alignment failed.\n' + str( e ) )
+ finally:
+ # clean up temp dir
+ if os.path.exists( tmp_index_dir ):
+ shutil.rmtree( tmp_index_dir )
+ if os.path.exists( tmp_dir ):
+ shutil.rmtree( tmp_dir )
+
if __name__=="__main__": __main__()
diff --git a/tools/sr_mapping/bwa_wrapper.xml b/tools/sr_mapping/bwa_wrapper.xml
index 39f82be598f..d90e34916de 100644
--- a/tools/sr_mapping/bwa_wrapper.xml
+++ b/tools/sr_mapping/bwa_wrapper.xml
@@ -27,23 +27,23 @@
--suppressHeader=$suppressHeader
-
-
-
-
-
+
+
+
+
+
-
-
-
-
-
-
+
+
+
+
+
+
-
-
+
+
-
+
@@ -64,11 +64,11 @@
-
-
+
+
-
+
@@ -78,7 +78,7 @@
-
+
@@ -94,6 +94,8 @@
bwa aln -t 4 phiX test-data/bwa_wrapper_in1.fastq > bwa_wrapper_out1.sai
bwa samse phiX bwa_wrapper_out1.sai test-data/bwa_wrapper_in1.fastq >> bwa_wrapper_out1.sam
phiX.fasta is the prefix for the reference
+ remove the comment lines (beginning with '@') from the resulting sam file
+ note that 'phiX' should be 'PHIX174' to match what's in the indexed file
-->
@@ -111,6 +113,7 @@
bwa aln -n 0.04 -o 1 -e -1 -d 16 -i 5 -k 2 -t 4 -M 3 -O 11 -E 4 -R -N phiX.fasta test-data/bwa_wrapper_in1.fastq > bwa_wrapper_out1.sai
bwa samse phiX.fasta bwa_wrapper_out1.sai test-data/bwa_wrapper_in1.fastq >> bwa_wrapper_out2.sam
phiX.fasta is the prefix for the reference
+ remove the comment lines (beginning with '@') from the resulting sam file
-->
@@ -118,7 +121,7 @@
-
+
@@ -143,6 +146,8 @@
bwa aln -n 0.04 -o 1 -e -1 -d 16 -i 5 -k 2 -t 4 -M 3 -O 11 -E 4 -R -N phiX.fasta test-data/bwa_wrapper_in3.fastq > bwa_wrapper_out3b.sai
bwa sampe -a 500 -o 100000 phiX.fasta bwa_wrapper_out3a.sai bwa_wrapper_out3b.sai test-data/bwa_wrapper_in2.fastq test-data/bwa_wrapper_in3.fastq >> bwa_wrapper_out3.sam
phiX.fasta is the prefix for the reference
+ remove the comment lines (beginning with '@') from the resulting sam file
+ note that 'phiX' should be 'PHIX174' to match what's in the indexed file
-->
@@ -150,8 +155,8 @@
-
-
+
+
@@ -171,12 +176,12 @@
-
-**What it does**
-
+
+**What it does**
+
BWA is a fast light-weighted tool that aligns relatively short sequences (queries) to a sequence database (large), such as the human reference genome. It is developed by Heng Li at the Sanger Insitute. Li H. and Durbin R. (2009) Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics, 25, 1754-60.
-This tool uses BWA version 0.5.3.
+This tool uses BWA version 0.5.5.
------
@@ -201,7 +206,7 @@ BWA accepts files in Sanger FASTQ format. Use the FASTQ Groomer to prepare your
The output is in SAM format, and has the following columns::
Column Description
- -------- --------------------------------------------------------
+ -------- --------------------------------------------------------
1 QNAME Query (pair) NAME
2 FLAG bitwise FLAG
3 RNAME Reference sequence NAME
@@ -250,53 +255,53 @@ This is an exhaustive list of BWA options:
For **aln**::
- -n NUM Maximum edit distance if the value is INT, or the fraction of missing
- alignments given 2% uniform base error rate if FLOAT. In the latter
+ -n NUM Maximum edit distance if the value is INT, or the fraction of missing
+ alignments given 2% uniform base error rate if FLOAT. In the latter
case, the maximum edit distance is automatically chosen for different
read lengths. [0.04]
-o INT Maximum number of gap opens [1]
- -e INT Maximum number of gap extensions, -1 for k-difference mode
+ -e INT Maximum number of gap extensions, -1 for k-difference mode
(disallowing long gaps) [-1]
-d INT Disallow a long deletion within INT bp towards the 3'-end [16]
-i INT Disallow an indel within INT bp towards the ends [5]
- -l INT Take the first INT subsequence as seed. If INT is larger than the
+ -l INT Take the first INT subsequence as seed. If INT is larger than the
query sequence, seeding will be disabled. For long reads, this option
is typically ranged from 25 to 35 for '-k 2'. [inf]
-k INT Maximum edit distance in the seed [2]
-t INT Number of threads (multi-threading mode) [1]
- -M INT Mismatch penalty. BWA will not search for suboptimal hits with a score
+ -M INT Mismatch penalty. BWA will not search for suboptimal hits with a score
lower than (bestScore-misMsc). [3]
-O INT Gap open penalty [11]
-E INT Gap extension penalty [4]
- -c Reverse query but not complement it, which is required for alignment
+ -c Reverse query but not complement it, which is required for alignment
in the color space.
- -R Proceed with suboptimal alignments even if the top hit is a repeat. By
- default, BWA only searches for suboptimal alignments if the top hit is
- unique. Using this option has no effect on accuracy for single-end
- reads. It is mainly designed for improving the alignment accuracy of
- paired-end reads. However, the pairing procedure will be slowed down,
+ -R Proceed with suboptimal alignments even if the top hit is a repeat. By
+ default, BWA only searches for suboptimal alignments if the top hit is
+ unique. Using this option has no effect on accuracy for single-end
+ reads. It is mainly designed for improving the alignment accuracy of
+ paired-end reads. However, the pairing procedure will be slowed down,
especially for very short reads (~32bp).
- -N Disable iterative search. All hits with no more than maxDiff
+ -N Disable iterative search. All hits with no more than maxDiff
differences will be found. This mode is much slower than the default.
For **samse**::
- -n INT Output up to INT top hits. Value -1 to disable outputting multiple
+ -n INT Output up to INT top hits. Value -1 to disable outputting multiple
hits. NOTE: Entering a value other than -1 will result in output that
- is not in SAM format, and therefore not usable further down the
- pipeline. Check the BWA documentation for details on the format of
+ is not in SAM format, and therefore not usable further down the
+ pipeline. Check the BWA documentation for details on the format of
the output. [-1]
For **sampe**::
- -a INT Maximum insert size for a read pair to be considered as being mapped
- properly. Since version 0.4.5, this option is only used when there
- are not enough good alignment to infer the distribution of insert
+ -a INT Maximum insert size for a read pair to be considered as being mapped
+ properly. Since version 0.4.5, this option is only used when there
+ are not enough good alignment to infer the distribution of insert
sizes. [500]
- -o INT Maximum occurrences of a read for pairing. A read with more
- occurrences will be treated as a single-end read. Reducing this
+ -o INT Maximum occurrences of a read for pairing. A read with more
+ occurrences will be treated as a single-end read. Reducing this
parameter helps faster pairing. [100000]
-
+