]> git.donarmstrong.com Git - samtools.git/commitdiff
make the script robust to the bugs in SOAP-2.1.7
authorHeng Li <lh3@live.co.uk>
Wed, 8 Apr 2009 09:57:08 +0000 (09:57 +0000)
committerHeng Li <lh3@live.co.uk>
Wed, 8 Apr 2009 09:57:08 +0000 (09:57 +0000)
misc/bowtie2sam.pl
misc/soap2sam.pl
misc/wgsim_eval.pl

index 4ead83991f81bab59c2a53485d2296e7b0ac6370..5dff88d4e96d25cba66e65cf4d3c5facaad22d21 100755 (executable)
@@ -12,7 +12,7 @@ exit;
 
 sub bowtie2sam {
   my %opts = ();
-  die("Usage: bowtie2sam.pl <aln.bowtie>\n") if (@ARGV == 0);
+  die("Usage: bowtie2sam.pl <aln.bowtie>\n") if (@ARGV == 0 && -t STDIN);
   # core loop
   my (@s, $last, @staging, $k, $best_s, $subbest_s, $best_k);
   $last = '';
index 5e721f912b783e2492e0428e30db5d3b56cab9ce..1628c0f2230b4cfbb45f3830e1e93bc2cb097e94 100755 (executable)
@@ -36,13 +36,14 @@ sub mating {
 sub soap2sam {
   my %opts = ();
   getopts("p", \%opts);
-  die("Usage: soap2sam.pl [-p] <aln.soap>\n") if (@ARGV == 0);
+  die("Usage: soap2sam.pl [-p] <aln.soap>\n") if (@ARGV == 0 && -t STDIN);
   my $is_paired = defined($opts{p});
   # core loop
   my @s1 = ();
   my @s2 = ();
   my ($s_last, $s_curr) = (\@s1, \@s2);
   while (<>) {
+       s/[\177-\377]|[\000-\010]|[\012-\040]//g;
        next if (&soap2sam_aux($_, $s_curr, $is_paired) < 0);
        if (@$s_last != 0 && $s_last->[0] eq $s_curr->[0]) {
          &mating($s_last, $s_curr);
@@ -60,8 +61,8 @@ sub soap2sam {
 sub soap2sam_aux {
   my ($line, $s, $is_paired) = @_;
   chomp($line);
-  my @t = split("\t", $line);
-  return -1 if (@t < 9);
+  my @t = split(/\s+/, $line);
+  return -1 if (@t < 9 || $line =~ /^\s/ || !$t[0]);
   @$s = ();
   # read name
   $s->[0] = $t[0];
@@ -71,7 +72,8 @@ sub soap2sam_aux {
   $s->[1] |= 1 | 1<<($t[4] eq 'a'? 6 : 7);
   $s->[1] |= 2 if ($is_paired);
   # read & quality
-  $s->[9] = $t[1]; $s->[10] = $t[2];
+  $s->[9] = $t[1];
+  $s->[10] = (length($t[2]) > length($t[1]))? substr($t[2], 0, length($t[1])) : $t[2];
   # cigar
   $s->[5] = length($s->[9]) . "M";
   # coor
index ecc53cca5163d8050ddb9159085d6b7f0ad3b9c3..99e2ac93b1c73fd37525bc963fcc4959a571cf9f 100755 (executable)
@@ -33,7 +33,6 @@ sub wgsim_eval {
        # skip unmapped reads
        next if (($t[1]&0x4) || $chr eq '*');
        # parse read name and check
-       ++$c0[$q];
        if ($t[0] =~ /^(\S+)_(\d+)_(\d+)_/) {
          if ($1 ne $chr) { # different chr
                $is_correct = 0;
@@ -57,8 +56,10 @@ sub wgsim_eval {
                }
          }
        } else {
-         die("[wgsim_eval] read '$t[0]' was not generated by wgsim?\n");
+         warn("[wgsim_eval] read '$t[0]' was not generated by wgsim?\n");
+         next;
        }
+       ++$c0[$q];
        ++$c1[$q] unless ($is_correct);
        print STDERR $line if (($flag&1) && !$is_correct && $q > 0);
   }