opensource/galaxy
3
0
Fork 0
mirror of https://github.com/galaxyproject/galaxy.git synced 2026-09-24 16:30:27 +08:00
Code Issues Packages Projects Releases Wiki Activity

Updated BWA wrapper, SAM to BAM, SAM Merge, and SAM Pileup so that the return code of the program is captured and if it fails, the stderr is printed; temp files are created locally; a message is printed if results file empty.

Browse Source
This commit is contained in:
Kelly Vincent
2010-02-23 10:40:19 -05:00
parent f47cd6a27f
commit a314b3716b
8 changed files with 418 additions and 219 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tools/samtools/sam_merge.py
+26 -9
View File
@@ -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__()
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__()
tools/samtools/sam_merge.xml
+42 -1
View File
@@ -19,7 +19,48 @@
<outputs>
<data name="output1" format="bam" />
</outputs>
<!-- bam files are binary and not sniffable so can't be uploaded without being corrupted, so no tests -->
<tests>
<!-- TODO: add ability to test framework to test without at least
one repeat element value
<test>
-->
<!--
Bam merge command:
samtools merge test-data/sam_merge_out1.bam test-data/sam_merge_in1.bam test-data/sam_merge_in2.bam
-->
<!--
<param name="input1" value="sam_merge_in1.bam" ftype="bam" />
<param name="input2" value="sam_merge_in2.bam" ftype="bam" />
<output name="output1" file="sam_merge_out1.bam" ftype="bam" />
</test>
-->
<test>
<!--
Bam merge command:
samtools merge test-data/sam_merge_out2.bam test-data/sam_merge_in1.bam test-data/sam_merge_in2.bam test-data/sam_merge_in3.bam
-->
<param name="input1" value="sam_merge_in1.bam" ftype="bam" />
<param name="input2" value="sam_merge_in2.bam" ftype="bam" />
<param name="input" value="sam_merge_in3.bam" ftype="bam" />
<output name="output1" file="sam_merge_out2.bam" ftype="bam" />
</test>
<!-- TODO: add ability to test code to be able to test with multiple
inputs (parameters with same value)
<test>
-->
<!--
Bam merge command:
samtools merge test-data/sam_merge_out3.bam test-data/sam_merge_in1.bam test-data/sam_merge_in2.bam test-data/sam_merge_in3.bam test-data/sam_merge_in4.bam
-->
<!--
<param name="input1" value="sam_merge_in1.bam" ftype="bam" />
<param name="input2" value="sam_merge_in2.bam" ftype="bam" />
<param name="input" value="sam_merge_in3.bam" ftype="bam" />
<param name="input" value="sam_merge_in4.bam" ftype="bam" />
<output name="output1" file="sam_merge_out3.bam" ftype="bam" />
</test>
-->
</tests>
<help>
**What it does**
tools/samtools/sam_pileup.py
+72 -58
View File
@@ -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__()
tools/samtools/sam_pileup.xml
+52 -17
View File
@@ -1,16 +1,16 @@
<tool id="sam_pileup" name="Generate pileup" version="1.0.0">
<description>from BAM dataset</description>
<command interpreter="python">
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 @@
<validator type="unspecified_build" />
<validator type="dataset_metadata_in_file" filename="sam_fa_indices.loc" metadata_name="dbkey" metadata_column="1" message="Sequences are not currently available for the specified build." line_startswith="index" />
</param>
</when>
</when>
<when value="history">
<param name="input1" type="data" format="bam" label="Select the BAM file to generate the pileup file for" />
<param name="ownFile" type="data" format="fasta" metadata_name="dbkey" label="Select a reference genome" />
</when>
</conditional>
</conditional>
<param name="lastCol" type="select" label="Whether or not to print the mapping quality as the last column" help="Makes the output easier to parse, but is space inefficient">
<option value="no">Do not print the mapping quality as the last column</option>
<option value="yes">Print the mapping quality as the last column</option>
@@ -72,10 +72,45 @@
<outputs>
<data format="tabular" name="output1" />
</outputs>
<!-- tests are not possible because bam is a non-sniffable binary format and cannot be uploaded without being corrupted -->
<tests>
<test>
<!--
Bam to pileup command:
samtools faidx chr_m.fasta
samtools pileup -M 60 -f chr_m.fasta test-data/sam_pileup_in1.bam > test-data/sam_pileup_out1.pileup
chr_m.fasta is the prefix of the index
-->
<param name="reference" value="history" />
<param name="input1" value="sam_pileup_in1.bam" ftype="bam" />
<param name="ownFile" value="chr_m.fasta" ftype="fasta" />
<param name="lastCol" value="no" />
<param name="indels" value="no" />
<param name="mapCap" value="60" />
<param name="consensus" value="no" />
<output name="output1" file="sam_pileup_out1.pileup" />
</test>
<test>
<!--
Bam to pileup command:
samtools pileup -M 60 -c -T 0.85 -N 2 -r 0.001 -I 40 -f chr_m.fasta test-data/sam_pileup_in1.bam > test-data/sam_pileup_out2.pileup
chr_m.fasta is the prefix of the index
-->
<param name="reference" value="indexed" />
<param name="input1" value="sam_pileup_in1.bam" ftype="bam" dbkey="chrM" />
<param name="lastCol" value="no" />
<param name="indels" value="no" />
<param name="mapCap" value="60" />
<param name="consensus" value="yes" />
<param name="theta" value="0.85" />
<param name="hapNum" value="2" />
<param name="fraction" value="0.001" />
<param name="phredProb" value="40" />
<output name="output1" file="sam_pileup_out2.pileup" />
</test>
</tests>
<help>
**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
tools/samtools/sam_to_bam.py
+62 -30
View File
@@ -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 <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()
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__()
tools/samtools/sam_to_bam.xml
+19
View File
@@ -32,11 +32,30 @@ sam_to_bam.py --input1=$source.input1 --dbkey=${input1.metadata.dbkey}
</outputs>
<tests>
<test>
<!--
Sam-to-Bam command:
samtools faidx chr_m.fasta
samtools view -bt chr_m.fasta.fai -o unsorted.bam test-data/3.sam
samtools sort unsorted.bam sam_to_bam_out1
chr_m.fasta is the reference file
-->
<param name="index_source" value="history" />
<param name="input1" value="3.sam" ftype="sam" />
<param name="ref_file" value="chr_m.fasta" ftype="fasta" />
<output name="output1" file="sam_to_bam_out1.bam" />
</test>
<test>
<!--
Sam-to-Bam command:
samtools view -bt chr_m.fasta.fai -o unsorted.bam test-data/3.sam
samtools sort unsorted.bam sam_to_bam_out2
chr_m.fasta is the reference file and the index chr_m.fasta.fai
should be in the same directory
-->
<param name="index_source" value="cached" />
<param name="input1" value="3.sam" ftype="sam" dbkey="chrM" />
<param name="output1" file="sam_to_bam_out2.bam" />
</test>
</tests>
<help>
tools/sr_mapping/bwa_wrapper.py
+93 -57
View File
@@ -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__()
tools/sr_mapping/bwa_wrapper.xml
+52 -47
View File
@@ -27,23 +27,23 @@
--suppressHeader=$suppressHeader
</command>
<inputs>
<conditional name="genomeSource">
<param name="refGenomeSource" type="select" label="Will you select a reference genome from your history or use a built-in index?">
<option value="indexed">Use a built-in index</option>
<option value="history">Use one from the history</option>
</param>
<conditional name="genomeSource">
<param name="refGenomeSource" type="select" label="Will you select a reference genome from your history or use a built-in index?">
<option value="indexed">Use a built-in index</option>
<option value="history">Use one from the history</option>
</param>
<when value="indexed">
<param name="indices" type="select" label="Select a reference genome">
<options from_file="bwa_index.loc">
<column name="value" index="1" />
<column name="name" index="0" />
</options>
</param>
<param name="indices" type="select" label="Select a reference genome">
<options from_file="bwa_index.loc">
<column name="value" index="1" />
<column name="name" index="0" />
</options>
</param>
</when>
<when value="history">
<param name="ownFile" type="data" format="fasta" metadata_name="dbkey" label="Select a reference from history" />
<when value="history">
<param name="ownFile" type="data" format="fasta" metadata_name="dbkey" label="Select a reference from history" />
</when>
</conditional>
</conditional>
<conditional name="paired">
<param name="sPaired" type="select" label="Is this library mate-paired?">
<option value="single">Single-end</option>
@@ -64,11 +64,11 @@
</param>
<when value="pre_set" />
<when value="full">
<param name="maxEditDist" type="integer" value="0" label="Maximum edit distance (-n)" help="Enter this value OR a fraction of missing alignments, not both" />
<param name="fracMissingAligns" type="float" value="0.04" label="Fraction of missing alignments given 2% uniform base error rate (-n)" help="Enter this value OR maximum edit distance, not both" />
<param name="maxEditDist" type="integer" value="0" label="Maximum edit distance (-n)" help="Enter this value OR a fraction of missing alignments, not both" />
<param name="fracMissingAligns" type="float" value="0.04" label="Fraction of missing alignments given 2% uniform base error rate (-n)" help="Enter this value OR maximum edit distance, not both" />
<param name="maxGapOpens" type="integer" value="1" label="Maximum number of gap opens (-o)" />
<param name="maxGapExtens" type="integer" value="-1" label="Maximum number of gap extensions (-e)" help="-1 for k-difference mode (disallowing long gaps)" />
<param name="disallowLongDel" type="integer" value="16" label="Disallow long deletion within [value] towards the 3'-end (-d)" />
<param name="disallowLongDel" type="integer" value="16" label="Disallow long deletion within [value] bp towards the 3'-end (-d)" />
<param name="disallowIndel" type="integer" value="5" label="Disallow insertion/deletion within [value] bp towards the end (-i)" />
<param name="seed" type="integer" value="-1" label="Number of first subsequences to take as seed (-l)" help="Enter -1 for infinity" />
<param name="maxEditDistSeed" type="integer" value="2" label="Maximum edit distance in the seed (-k)" />
@@ -78,7 +78,7 @@
<param name="suboptAlign" type="boolean" truevalue="true" falsevalue="false" checked="no" label="Proceed with suboptimal alignments even if the top hit is a repeat" help="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) (-R)" />
<param name="noIterSearch" type="boolean" truevalue="true" falsevalue="false" checked="no" label="Disable iterative search" help="All hits with no more than maxDiff differences will be found. This mode is much slower than the default (-N)" />
<param name="outputTopN" type="integer" value="-1" label="Output top [value] hits" help="For single-end reads only. Enter -1 to disable outputting multiple hits. NOTE: If you put in a positive value here, your output will NOT be in SAM format (-n)" />
<param name="maxInsertSize" type="integer" value="500" label="Maximum insert size for a read pair to be considered as being mapped properly" help="For paired-end reads only. Only used when there are not enough good alignment to infer the distribution of insert sizes (-a)" />
<param name="maxInsertSize" type="integer" value="500" label="Maximum insert size for a read pair to be considered as being mapped properly" help="For paired-end reads only. Only used when there are not enough good alignments to infer the distribution of insert sizes (-a)" />
<param name="maxOccurPairing" type="integer" value="100000" label="Maximum occurrences of a read for pairing" help="For paired-end reads only. A read with more occurrences will be treated as a single-end read. Reducing this parameter helps faster pairing (-o)" />
</when>
</conditional>
@@ -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
-->
<param name="refGenomeSource" value="indexed" />
<param name="indices" value="phiX" />
@@ -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
-->
<param name="refGenomeSource" value="history" />
<param name="ownFile" value="phiX.fasta" />
@@ -118,7 +121,7 @@
<param name="input1" value="bwa_wrapper_in1.fastq" ftype="fastqsanger" />
<param name="source_select" value="full" />
<param name="maxEditDist" value="0" />
<param name="fracMissingAligns" value="0.04" />
<param name="fracMissingAligns" value="0.04" />
<param name="maxGapOpens" value="1" />
<param name="maxGapExtens" value="-1" />
<param name="disallowLongDel" value="16" />
@@ -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
-->
<param name="refGenomeSource" value="indexed" />
<param name="indices" value="phiX" />
@@ -150,8 +155,8 @@
<param name="input1" value="bwa_wrapper_in2.fastq" ftype="fastqsanger" />
<param name="input2" value="bwa_wrapper_in3.fastq" ftype="fastqsanger" />
<param name="source_select" value="full" />
<param name="maxEditDist" value="0" />
<param name="fracMissingAligns" value="0.04" />
<param name="maxEditDist" value="0" />
<param name="fracMissingAligns" value="0.04" />
<param name="maxGapOpens" value="1" />
<param name="maxGapExtens" value="-1" />
<param name="disallowLongDel" value="16" />
@@ -171,12 +176,12 @@
</test>
</tests>
<help>
**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]
</help>
<code file="bwa_wrapper_code.py" />
</tool>