From 47630c9ed6e32ea2bb25abaaca7582f995e65752 Mon Sep 17 00:00:00 2001
From: Andrea Sbardellati <andrea.sbardellati@crs4.it>
Date: Fri, 18 Jan 2013 17:11:04 +0100
Subject: [PATCH] Add option for writing not aligned reads to file

Reviewed-by: Nicola Soranzo <soranzo@crs4.it>
---
 tools/sr_mapping/bowtie2_wrapper.py  | 33 ++++++++++++++++++++++++++--
 tools/sr_mapping/bowtie2_wrapper.xml | 32 +++++++++++++++++++++++++--
 2 files changed, 61 insertions(+), 4 deletions(-)

diff --git a/tools/sr_mapping/bowtie2_wrapper.py b/tools/sr_mapping/bowtie2_wrapper.py
index 496f4786a7a..807faaf23f9 100644
--- a/tools/sr_mapping/bowtie2_wrapper.py
+++ b/tools/sr_mapping/bowtie2_wrapper.py
@@ -31,10 +31,29 @@ def __main__():
     parser.add_option( '', '--rgpl', dest='rgpl', help='Platform/technology used to produce the reads' )
     parser.add_option( '', '--rgsm', dest='rgsm', help='Sample' )
 
+    parser.add_option( '', '--output_unaligned_reads', dest='output_unaligned_reads', help='File name for unaligned reads (single-end)' )
+    parser.add_option( '', '--output_unaligned_reads_l', dest='output_unaligned_reads_l', help='File name for unaligned reads (left, paired-end)' )
+    parser.add_option( '', '--output_unaligned_reads_r', dest='output_unaligned_reads_r', help='File name for unaligned reads (right, paired-end)' )
+
     (options, args) = parser.parse_args()
     
+    tmp_unaligned_file_name = None
     # Creat bowtie index if necessary.
     tmp_index_dir = tempfile.mkdtemp()
+    
+    if options.single_paired == 'paired':
+        if options.output_unaligned_reads_l and options.output_unaligned_reads_r:
+            tmp_unaligned_file = tempfile.NamedTemporaryFile( dir=tmp_index_dir, suffix='.fastq' )
+            tmp_unaligned_file_name = tmp_unaligned_file.name
+            tmp_unaligned_file.close()
+            output_unaligned_reads = '--un-conc %s' % tmp_unaligned_file_name
+        else:
+            output_unaligned_reads = ''
+    elif options.output_unaligned_reads:
+        output_unaligned_reads = '--un %s' % options.output_unaligned_reads
+    else:
+        output_unaligned_reads = ''
+    
     if options.own_file:
         index_path = os.path.join( tmp_index_dir, '.'.join( os.path.split( options.own_file )[1].split( '.' )[:-1] ) )
         try:
@@ -71,7 +90,7 @@ def __main__():
         index_path = options.index_path
 
     # Build bowtie command; use view and sort to create sorted bam.
-    cmd = 'bowtie2 %s -x %s %s | samtools view -Su - | samtools sort -o - - > %s'
+    cmd = 'bowtie2 %s -x %s %s %s | samtools view -Su - | samtools sort -o - - > %s'
     
     # Set up reads.
     if options.single_paired == 'paired':
@@ -104,7 +123,7 @@ def __main__():
         opts += ' --rg %s:%s' % ( 'SM', options.rgsm )
         
     # Final command:
-    cmd = cmd % ( opts, index_path, reads, options.output )
+    cmd = cmd % ( opts, index_path, reads, output_unaligned_reads, options.output )
     print cmd
 
     # Run
@@ -136,6 +155,16 @@ def __main__():
     except Exception, e:
         stop_err( 'Error in bowtie2:\n' + str( e ) ) 
 
+    # get unaligned reads output files in place if appropriate
+    if options.single_paired == 'paired' and tmp_unaligned_file_name and options.output_unaligned_reads_l and options.output_unaligned_reads_r:
+        try:
+            left = tmp_unaligned_file_name.replace( '.fastq', '.1.fastq' )
+            right = tmp_unaligned_file_name.replace( '.fastq', '.2.fastq' )
+            shutil.move( left, options.output_unaligned_reads_l )
+            shutil.move( right, options.output_unaligned_reads_r )
+        except Exception, e:
+            sys.stdout.write( 'Error producing the unaligned output files.\n' )
+
     # Clean up temp dirs
     if os.path.exists( tmp_index_dir ):
         shutil.rmtree( tmp_index_dir )
diff --git a/tools/sr_mapping/bowtie2_wrapper.xml b/tools/sr_mapping/bowtie2_wrapper.xml
index 8d444ce2c5c..6786c3eb761 100644
--- a/tools/sr_mapping/bowtie2_wrapper.xml
+++ b/tools/sr_mapping/bowtie2_wrapper.xml
@@ -1,6 +1,6 @@
 <tool id="bowtie2" name="Bowtie2" version="0.1">
     <!-- Wrapper compatible with Bowtie version 2.0.0 -->
-    <description>is a short-read mapper</description>
+    <description>is a short-read aligner</description>
     <version_command>bowtie2 --version</version_command>
     <requirements>
         <requirement type="package">bowtie2</requirement>
@@ -15,6 +15,17 @@
         ## Outputs.
         --output=$output
 
+        #if str( $singlePaired.sPaired ) == "single"
+            #if $output_unaligned_reads_l
+                --output_unaligned_reads=$output_unaligned_reads_l
+            #end if
+        #else
+            #if $output_unaligned_reads_l and $output_unaligned_reads_r
+                --output_unaligned_reads_l=$output_unaligned_reads_l
+                --output_unaligned_reads_r=$output_unaligned_reads_r
+            #end if
+        #end if
+        
         ## Handle reference file.
         #if $refGenomeSource.genomeSource == "history":
             --own-file=$refGenomeSource.ownFile
@@ -71,6 +82,7 @@
                 <param name="maxInsert" type="integer" value="250" label="Maximum insert size for valid paired-end alignments" />
             </when>
         </conditional>
+        <param name="unalignedFile" type="boolean" truevalue="true" falsevalue="false" checked="False" label="Write all reads that could not be aligned to a file (uses --un for single-end and  --un-conc for paired-ends)" />
         <conditional name="refGenomeSource">
           <param name="genomeSource" type="select" label="Will you select a reference genome from your history or use a built-in index?" help="Built-ins were indexed using default options">
             <option value="indexed">Use a built-in index</option>
@@ -133,7 +145,7 @@
     </inputs>
 
     <outputs>
-        <data format="bam" name="output" label="${tool.name} on ${on_string}: mapped reads">
+        <data format="bam" name="output" label="${tool.name} on ${on_string}: aligned reads">
           <actions>
             <conditional name="refGenomeSource.genomeSource">
               <when value="indexed">
@@ -152,6 +164,22 @@
             </conditional>
           </actions>
         </data>
+        <data format="fastqsanger" name="output_unaligned_reads_l" label="${tool.name} on ${on_string}: unaligned reads (L)" >
+            <filter>unalignedFile is True</filter>
+            <actions>
+                <action type="format">
+                    <option type="from_param" name="singlePaired.input1" param_attribute="ext" />
+                </action>
+            </actions>
+        </data>
+        <data format="fastqsanger" name="output_unaligned_reads_r" label="${tool.name} on ${on_string}: unaligned reads (R)">
+            <filter>singlePaired['sPaired'] == "paired" and unalignedFile is True</filter>
+            <actions>
+                <action type="format">
+                    <option type="from_param" name="singlePaired.input1" param_attribute="ext" />
+                </action>
+            </actions>
+        </data>
     </outputs>
 
     <tests>