c_cpp/lib/htslib/test/compare_sam.pl

   1 #!/usr/bin/perl -w
   2 #
   3 #    Copyright (C) 2013 Genome Research Ltd.
   4 #
   5 #    Author: James Bonfield <jkb@sanger.ac.uk>
   6 #
   7 # Permission is hereby granted, free of charge, to any person obtaining a copy
   8 # of this software and associated documentation files (the "Software"), to deal
   9 # in the Software without restriction, including without limitation the rights
  10 # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  11 # copies of the Software, and to permit persons to whom the Software is
  12 # furnished to do so, subject to the following conditions:
  13 #
  14 # The above copyright notice and this permission notice shall be included in
  15 # all copies or substantial portions of the Software.
  16 #
  17 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  18 # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  19 # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  20 # THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  21 # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  22 # FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  23 # DEALINGS IN THE SOFTWARE.
  24
  25 # Compares two SAM files to report differences.
  26 # Optionally can skip header or ignore specific types of diff.
  27
  28 use strict;
  29 use Getopt::Long;
  30
  31 my %opts;
  32 GetOptions(\%opts, 'noqual', 'noaux', 'notemplate', 'unknownrg', 'nomd', 'template-1', 'noflag', 'Baux');
  33
  34 my ($fn1, $fn2) = @ARGV;
  35 open(my $fd1, "<", $fn1) || die $!;
  36 open(my $fd2, "<", $fn2) || die $!;
  37
  38 # Headers
  39 my ($c1,$c2)=(1,1);
  40 my (@hd1, @hd2, $ln1, $ln2);
  41 while (<$fd1>) {
  42     if (/^@/) {
  43         push(@hd1, $_);
  44     } else {
  45         $ln1 = $_;
  46         last;
  47     }
  48     $c1++;
  49 }
  50
  51 while (<$fd2>) {
  52     if (/^@/) {
  53         push(@hd2, $_);
  54     } else {
  55         $ln2 = $_;
  56         last;
  57     }
  58     $c2++;
  59 }
  60
  61 # FIXME: to do
  62 #print "@hd1\n";
  63 #print "@hd2\n";
  64
  65 # Compare lines
  66 while ($ln1 && $ln2) {
  67     chomp($ln1);
  68     chomp($ln2);
  69
  70     # Java CRAM adds RG:Z:UNKNOWN when the read-group is absent
  71     if (exists $opts{unknownrg}) {
  72         $ln1 =~ s/\tRG:Z:UNKNOWN//;
  73         $ln2 =~ s/\tRG:Z:UNKNOWN//;
  74     }
  75
  76     if (exists $opts{nomd}) {
  77         $ln1 =~ s/\tMD:Z:[A-Z0-9^]*//;
  78         $ln2 =~ s/\tMD:Z:[A-Z0-9^]*//;
  79         $ln1 =~ s/\tNM:i:\d+//;
  80         $ln2 =~ s/\tNM:i:\d+//;
  81     }
  82
  83     my @ln1 = split("\t", $ln1);
  84     my @ln2 = split("\t", $ln2);
  85
  86     # Fix BWA bug: unmapped data should have no alignments
  87     if ($ln1[1] & 4) { $ln1[4] = 0; $ln1[5] = "*"; }
  88     if ($ln2[1] & 4) { $ln2[4] = 0; $ln2[5] = "*"; }
  89
  90     # Canonicalise floating point numbers
  91     map {s/^(..):f:(.*)/{"$1:f:".($2+0)}/e} @ln1[11..$#ln1];
  92     map {s/^(..):f:(.*)/{"$1:f:".($2+0)}/e} @ln2[11..$#ln2];
  93
  94
  95     if (exists $opts{Baux}) {
  96         # Turn ??:H:<hex> into ??:B:c,<vals> so we can compare
  97         # Cramtools.jar vs htslib encodings.  Probably doable with (un)pack
  98         map {s/^(..):H:(.*)/{join(",", "$1:B:C", map {hex $_} $2=~m:..:g)}/e} @ln1[11..$#ln1];
  99         map {s/^(..):H:(.*)/{join(",", "$1:B:C", map {hex $_} $2=~m:..:g)}/e} @ln2[11..$#ln2];
 100
 101         # Canonicalise ??:B:? data series to be unsigned
 102         map {s/^(..):B:c(,?)(.*)/{"$1:B:C$2".join(",",map {($_+256)&255} split(",",$3))}/e} @ln1[11..$#ln1];
 103         map {s/^(..):B:c(,?)(.*)/{"$1:B:C$2".join(",",map {($_+256)&255} split(",",$3))}/e} @ln2[11..$#ln2];
 104
 105         map {s/^(..):B:s(,?)(.*)/{"$1:B:S$2".join(",",map {($_+65536)&65535} split(",",$3))}/e} @ln1[11..$#ln1];
 106         map {s/^(..):B:s(,?)(.*)/{"$1:B:S$2".join(",",map {($_+65536)&65535} split(",",$3))}/e} @ln2[11..$#ln2];
 107
 108         map {s/^(..):B:i(,?)(.*)/{"$1:B:I$2".join(",",map {$_<0? ($_+4294967296) : $_} split(",",$3))}/e} @ln1[11..$#ln1];
 109         map {s/^(..):B:i(,?)(.*)/{"$1:B:I$2".join(",",map {$_<0? ($_+4294967296) : $_} split(",",$3))}/e} @ln2[11..$#ln2];
 110     }
 111
 112     # Rationalise order of auxiliary fields
 113     if (exists $opts{noaux}) {
 114         @ln1 = @ln1[0..10];
 115         @ln2 = @ln2[0..10];
 116     } else {
 117         #my @a=@ln1[11..$#ln1];print "<<<@a>>>\n";
 118         @ln1[11..$#ln1] = sort @ln1[11..$#ln1];
 119         @ln2[11..$#ln2] = sort @ln2[11..$#ln2];
 120     }
 121
 122     if (exists $opts{noqual}) {
 123         $ln1[10] = "*";
 124         $ln2[10] = "*";
 125     }
 126
 127     if (exists $opts{notemplate}) {
 128         @ln1[6..8] = qw/* 0 0/;
 129         @ln2[6..8] = qw/* 0 0/;
 130     }
 131
 132     if (exists $opts{noflag}) {
 133         $ln1[1] = 0; $ln2[1] = 0;
 134     }
 135
 136     if (exists $opts{'template-1'}) {
 137         if (abs($ln1[8] - $ln2[8]) == 1) {
 138             $ln1[8] = $ln2[8];
 139         }
 140     }
 141
 142     # Cram doesn't uppercase the reference
 143     $ln1[9] = uc($ln1[9]);
 144     $ln2[9] = uc($ln2[9]);
 145
 146     # Cram will populate a sequence string that starts as "*"
 147     $ln2[9] = "*" if ($ln1[9] eq "*");
 148
 149     # Fix 0<op> cigar fields
 150     $ln1[5] =~ s/(\D|^)0\D/$1/g;
 151     $ln1[5] =~ s/^$/*/g;
 152     $ln2[5] =~ s/(\D|^)0\D/$1/g;
 153     $ln2[5] =~ s/^$/*/g;
 154
 155     # Fix 10M10M cigar to 20M
 156     $ln1[5] =~ s/(\d+)(\D)(\d+)(\2)/$1+$3.$2/e;
 157     $ln2[5] =~ s/(\d+)(\D)(\d+)(\2)/$1+$3.$2/e;
 158
 159     if ("@ln1" ne "@ln2") {
 160         print "Diff at lines $fn1:$c1, $fn2:$c2\n";
 161         my @s1 = split("","@ln1");
 162         my @s2 = split("","@ln2");
 163         my $ptr = "";
 164         for (my $i=0; $i < $#s1; $i++) {
 165             if ($s1[$i] eq $s2[$i]) {
 166                 $ptr .= "-";
 167             } else {
 168                 last;
 169             }
 170         }
 171         print "1\t@ln1\n2\t@ln2\n\t$ptr^\n\n";
 172         exit(1);
 173     }
 174
 175     $ln1 = <$fd1>;
 176     $ln2 = <$fd2>;
 177
 178     $c1++; $c2++;
 179 }
 180
 181 if (defined($ln1)) {
 182     print "EOF on $fn1\n";
 183     exit(1);
 184 }
 185
 186 if (defined($ln2)) {
 187     print "EOF on $fn2\n";
 188     exit(1);
 189 }
 190
 191 close($fd1);
 192 close($fd2);
 193
 194 exit(0);