my $flt = 0;
# parse annotations
my ($dp, $mq, $dp_alt) = (-1, -1, -1);
- if ($t[7] =~ /DP=(\d+)/i) {
- $dp = $1;
- } elsif ($t[7] =~ /DP4=(\d+),(\d+),(\d+),(\d+)/i) {
+ if ($t[7] =~ /DP4=(\d+),(\d+),(\d+),(\d+)/i) {
$dp = $1 + $2 + $3 + $4;
$dp_alt = $3 + $4;
}
+ if ($t[7] =~ /DP=(\d+)/i) {
+ $dp = $1;
+ }
$mq = $1 if ($t[7] =~ /MQ=(\d+)/i);
# the depth and mapQ filter
if ($dp >= 0) {
my @s = split(',', $t[4]);
for my $x (@s) {
my $l = length($x) - length($t[3]) + 5000;
+ if ($x =~ /^-/) {
+ $l = -(length($x) - 1) + 5000;
+ } elsif ($x =~ /^\+/) {
+ $l = length($x) - 1 + 5000;
+ }
$c0[$l] += 1 / @s;
}
}
for (my $i = 0; $i < 10000; ++$i) {
next if ($c0[$i] == 0);
- printf("%d\t%.2f\n", ($i-5000), $c0[$i]);
+ $c1[0] += $c0[$i];
+ $c1[1] += $c0[$i] if (($i-5000)%3 == 0);
+ printf("C\t%d\t%.2f\n", ($i-5000), $c0[$i]);
}
+ printf("3\t%d\t%d\t%.3f\n", $c1[0], $c1[1], $c1[1]/$c1[0]);
}
sub ucscsnp2vcf {