Added header removal option to BWA wrapper tool, and cleaned up its test files. Also added the bwa_solid2fastq_modified.pl script to the next gen conversion directory for the solid-to-fastq tool.

tools/next_gen_conversion/bwa_solid2fastq_modified.pl

@@ -0,0 +1,112 @@
+#!/usr/bin/perl -w
+
+# Author: lh3
+# Note: Ideally, this script should be written in C. It is a bit slow at present.
+
+use strict;
+use warnings;
+use Getopt::Std;
+
+my %opts;
+my $version = '0.1.2';
+my $usage = qq{
+Usage: solid2fastq.pl <paired> <outfile1> <outfile2> <outfile3> <F3.csfasta> <F3.qual> <R3.csfasta> <R3.qual> 
+
+Note: <in.title> is the string showed in the `# Title:' line of a
+      ".csfasta" read file. Then <in.title>F3.csfasta is read sequence
+      file and <in.title>F3_QV.qual is the quality file. If
+      <in.title>R3.csfasta is present, this script assumes reads are
+      paired; otherwise reads will be regarded as single-end.
+
+      The read name will be <out.prefix>:panel_x_y/[12] with `1' for R3
+      tag and `2' for F3. Usually you may want to use short <out.prefix>
+      to save diskspace. Long <out.prefix> also causes troubles to maq.
+
+};
+
+getopts('', \%opts);
+die($usage) if (@ARGV != 8);
+my ($is_paired,$outfile1,$outfile2,$outfile3,$f3reads,$f3qual,$r3reads,$r3qual) = @ARGV;
+my (@fhr, @fhw);
+my $fn = '';
+my @fn_suff = ($f3reads,$f3qual,$r3reads,$r3qual);
+#my @fn_suff = ('F3.csfasta', 'F3_QV.qual', 'R3.csfasta', 'R3_QV.qual');
+#my $is_paired = (-f "$title$fn_suff[2]" || -f "$title$fn_suff[2].gz")? 1 : 0;
+if ($is_paired eq "yes") { # paired end
+  for (0 .. 3) {
+	$fn = $fn_suff[$_];
+	$fn = "gzip -dc $fn.gz |" if (!-f $fn && -f "$fn.gz");
+	open($fhr[$_], $fn) || die("** Fail to open '$fn'.\n");
+  }
+  open($fhw[0], "|gzip >$outfile2") || die;
+  open($fhw[1], "|gzip >$outfile1") || die;
+  open($fhw[2], "|gzip >$outfile3") || die;
+  my (@df, @dr);
+  @df = &read1(1); @dr = &read1(2);
+  while (@df && @dr) {
+	if ($df[0] eq $dr[0]) { # mate pair
+	  print {$fhw[0]} $df[1]; print {$fhw[1]} $dr[1];
+	  @df = &read1(1); @dr = &read1(2);
+	} else {
+	  if ($df[0] le $dr[0]) {
+		print {$fhw[2]} $df[1];
+		@df = &read1(1);
+	  } else {
+		print {$fhw[2]} $dr[1];
+		@dr = &read1(2);
+	  }
+	}
+  }
+  if (@df) {
+	print {$fhw[2]} $df[1];
+	while (@df = &read1(1, $fhr[0], $fhr[1])) {
+	  print {$fhw[2]} $df[1];
+	}
+  }
+  if (@dr) {
+	print {$fhw[2]} $dr[1];
+	while (@dr = &read1(2, $fhr[2], $fhr[3])) {
+	  print {$fhw[2]} $dr[1];
+	}
+  }
+  close($fhr[$_]) for (0 .. $#fhr);
+  close($fhw[$_]) for (0 .. $#fhw);
+} else { # single end
+  for (0 .. 1) {
+	my $fn = "$fn_suff[$_]";
+	$fn = "gzip -dc $fn.gz |" if (!-f $fn && -f "$fn.gz");
+	open($fhr[$_], $fn) || die("** Fail to open '$fn'.\n");
+  }
+  open($fhw[2], "|gzip >$outfile1") || die;
+  my @df;
+  while (@df = &read1(1, $fhr[0], $fhr[1])) {
+	print {$fhw[2]} $df[1];
+  }
+  close($fhr[$_]) for (0 .. $#fhr);
+  close($fhw[2]);
+}
+
+sub read1 {
+  my $i = shift(@_);
+  my $j = ($i-1)<<1;
+  my ($key, $seq);
+  my ($fhs, $fhq) = ($fhr[$j], $fhr[$j|1]);
+  while (<$fhs>) {
+	my $t = <$fhq>;
+	if (/^>(\d+)_(\d+)_(\d+)_[FR]3/) {
+	  $key = sprintf("%.4d_%.4d_%.4d", $1, $2, $3); # this line could be improved on 64-bit machines
+	  #print  $key;
+	  die(qq/** unmatched read name: '$_' != '$_'\n/) unless ($_ eq $t);
+	  my $name = "$1_$2_$3/$i";
+	  $_ = substr(<$fhs>, 2);
+	  tr/0123./ACGTN/;
+	  my $s = $_;
+	  $_ = <$fhq>;
+	  s/^(\d+)\s*//;
+	  s/(\d+)\s*/chr($1+33)/eg;
+	  $seq = qq/\@$name\n$s+\n$_\n/;
+	  last;
+	}
+  }
+  return defined($seq)? ($key, $seq) : ();
+}

tools/sr_mapping/bwa_wrapper.py

@@ -41,6 +41,7 @@ def __main__():
     parser.add_option('', '--maxInsertSize', dest='maxInsertSize', help='Maximum insert size for a read pair to be considered mapped good')
     parser.add_option('', '--maxOccurPairing', dest='maxOccurPairing', help='Maximum occurrences of a read for pairings')
     parser.add_option('', '--dbkey', dest='dbkey', help='')
+    parser.add_option('', '--suppressHeader', dest='suppressHeader', help='Suppress header')
     (options, args) = parser.parse_args()
        
     # index if necessary
@@ -119,5 +120,29 @@ def __main__():
     # clean up temp files
     tmp_align_out.close()
     tmp_align_out2.close()
+    # remove header if necessary
+    if options.suppressHeader == 'true':
+        tmp_out = tempfile.NamedTemporaryFile()
+        cmd4 = 'cp %s %s' % (options.output, tmp_out.name)
+        try:
+            os.system(cmd4)
+        except Exception, erf:
+            stop_err("Error copying output file before removing headers\n" + str(erf))
+        output = file(tmp_out.name, 'r')
+        fout = file(options.output, 'w')
+        header = True
+        line = output.readline()
+        while line.strip() != '':
+            if header:
+                if line.startswith('@HD') or line.startswith('@SQ') or line.startswith('@RG') or line.startswith('@PG') or line.startswith('@CO'):
+                    pass
+                else:
+                    header = False
+                    fout.write(line)
+            else:
+                fout.write(line)
+            line = output.readline()
+        fout.close()
+        tmp_out.close()
     
 if __name__=="__main__": __main__()

tools/sr_mapping/bwa_wrapper.xml

@@ -61,6 +61,7 @@
     #else:
      --dbkey="None"
     #end if
+    --suppressHeader=$suppressHeader
   </command>
   <inputs>
     <conditional name="solidOrSolexa">
@@ -151,6 +152,7 @@
         <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" />
       </when>
     </conditional>
+    <param name="suppressHeader" type="boolean" truevalue="true" falsevalue="false" checked="true" label="Suppress the header in the output SAM file" help="BWA produces SAM with several lines of header information" />
   </inputs>
   <outputs>
     <data format="sam" name="output" />
@@ -163,6 +165,7 @@
       <param name="sPaired" value="single" />
       <param name="input1" value="bwa_wrapper_in1.fastq" />
       <param name="source_select" value="pre_set" />
+      <param name="suppressHeader" value="true" />
       <output name="output" file="bwa_wrapper_out0.sam" ftype="sam" />
     </test>
     <test>
@@ -172,6 +175,7 @@
       <param name="sPaired" value="single" />
       <param name="input1" value="bwa_wrapper_in0.fastq" />
       <param name="source_select" value="pre_set" />
+      <param name="suppressHeader" value="true" />
       <output name="output" file="bwa_wrapper_out0b.sam" ftype="sam" />
     </test>
     <test>
@@ -198,6 +202,7 @@
       <param name="outputTopN" value="-1" />
       <param name="maxInsertSize" value="500" />
       <param name="maxOccurPairing" value="100000" />
+      <param name="suppressHeader" value="true" />
       <output name="output" file="bwa_wrapper_out1.sam" ftype="sam" />
     </test>
     <test>
@@ -225,6 +230,7 @@
       <param name="outputTopN" value="-1" />
       <param name="maxInsertSize" value="500" />
       <param name="maxOccurPairing" value="100000" />
+      <param name="suppressHeader" value="true" />
       <output name="output" file="bwa_wrapper_out2.sam" ftype="sam" />
     </test> 
     <test>
@@ -251,6 +257,7 @@
       <param name="outputTopN" value="-1" />
       <param name="maxInsertSize" value="500" />
       <param name="maxOccurPairing" value="100000" />
+      <param name="suppressHeader" value="true" />
       <output name="output" file="bwa_wrapper_out3.sam" ftype="sam" />
     </test> 
     <test>
@@ -278,8 +285,9 @@
       <param name="outputTopN" value="-1" />
       <param name="maxInsertSize" value="500" />
       <param name="maxOccurPairing" value="100000" />
+      <param name="suppressHeader" value="true" />
       <output name="output" file="bwa_wrapper_out4.sam" ftype="sam" />
-    </test>
+    </test> 
   </tests> 
   <help>