$_ = $t[5]; s/(\d+)[MDN]/$rght+=$1,'x'/eg;
--$rght;
# correct for soft clipping
+ my ($left0, $rght0) = ($left, $rght);
$left -= $1 if (/^(\d+)[SH]/);
$rght += $1 if (/(\d+)[SH]$/);
+ $left0 -= $1 if (/(\d+)[SH]$/);
+ $rght0 += $1 if (/^(\d+)[SH]/);
# skip unmapped reads
next if (($t[1]&0x4) || $chr eq '*');
# parse read name and check
} else {
if ($flag & 2) {
if (($t[1]&0x40) && !($t[1]&0x10)) { # F3, forward
- $is_correct = 0 if (abs($2 - $left) > $gap);
+ $is_correct = 0 if (abs($2 - $left) > $gap && abs($2 - $left0) > $gap);
} elsif (($t[1]&0x40) && ($t[1]&0x10)) { # F3, reverse
- $is_correct = 0 if (abs($3 - $rght) > $gap);
+ $is_correct = 0 if (abs($3 - $rght) > $gap && abs($3 - $rght0) > $gap);
} elsif (($t[1]&0x80) && !($t[1]&0x10)) { # R3, forward
- $is_correct = 0 if (abs($3 - $left) > $gap);
+ $is_correct = 0 if (abs($3 - $left) > $gap && abs($3 - $left0) > $gap);
} else { # R3, reverse
- $is_correct = 0 if (abs($2 - $rght) > $gap);
+ $is_correct = 0 if (abs($2 - $rght) > $gap && abs($3 - $rght0) > $gap);
}
} else {
if ($t[1] & 0x10) { # reverse
- $is_correct = 0 if (abs($3 - $rght) > $gap); # in case of indels that are close to the end of a reads
+ $is_correct = 0 if (abs($3 - $rght) > $gap && abs($3 - $rght0) > $gap); # in case of indels that are close to the end of a reads
} else {
- $is_correct = 0 if (abs($2 - $left) > $gap);
+ $is_correct = 0 if (abs($2 - $left) > $gap && abs($2 - $left0) > $gap);
}
}
}