annotate ezBAMQC/src/htslib/test/compare_sam.pl @ 18:494b5cd02238

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