diff --git a/tools/next_gen_conversion/bwa_solid2fastq_modified.pl b/tools/next_gen_conversion/bwa_solid2fastq_modified.pl new file mode 100755 index 00000000000..d38d38732e1 --- /dev/null +++ b/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 + +Note: is the string showed in the `# Title:' line of a + ".csfasta" read file. Then F3.csfasta is read sequence + file and F3_QV.qual is the quality file. If + R3.csfasta is present, this script assumes reads are + paired; otherwise reads will be regarded as single-end. + + The read name will be :panel_x_y/[12] with `1' for R3 + tag and `2' for F3. Usually you may want to use short + to save diskspace. Long 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) : (); +} diff --git a/tools/sr_mapping/bwa_wrapper.py b/tools/sr_mapping/bwa_wrapper.py index 82652e2e2d6..1bece228a2d 100644 --- a/tools/sr_mapping/bwa_wrapper.py +++ b/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__() diff --git a/tools/sr_mapping/bwa_wrapper.xml b/tools/sr_mapping/bwa_wrapper.xml index 10d806c28f2..021f367145a 100644 --- a/tools/sr_mapping/bwa_wrapper.xml +++ b/tools/sr_mapping/bwa_wrapper.xml @@ -61,6 +61,7 @@ #else: --dbkey="None" #end if + --suppressHeader=$suppressHeader @@ -151,6 +152,7 @@ + @@ -163,6 +165,7 @@ + @@ -172,6 +175,7 @@ + @@ -198,6 +202,7 @@ + @@ -225,6 +230,7 @@ + @@ -251,6 +257,7 @@ + @@ -278,8 +285,9 @@ + - +