1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82
|
#!/usr/bin/env perl
# Show SNPs around an INDEL
# for line in the vcf
# stuff the line into a queue
# when you reach an indel
# record the position
# pop lines from the back of the queue until we are at the current position
#
my @lines;
my $prox = $ARGV[0];
my $lastchrom = "";
my $indelpos = 0;
while (<STDIN>) {
if ($_ =~ /^#/) {
print $_;
next;
}
$_ =~ /^(.+?)\t(.+?)\t(.+?)\t(.+?)\t(.+?)\t/;
my $chrom = $1;
my $pos = $2;
my $tag = $3;
my $ref = $4;
my $alt = $5;
#print "chrom: $chrom, pos: $pos, ref: $ref, alt: $alt\n";
# if new chrom, print out everything from last one
if ($lastchrom == "") {
$lastchrom = $chrom;
}
if ($chrom != $lastchrom) {
while ($lines) {
print pop(@lines);
}
}
unshift(@lines, $_);
my $diff = length($ref) - length($alt);
if ($diff != 0) {
# insertion
if ($indelpos == 0) {
$indelpos = $pos;
}
$nextindelpos = $pos;
#print "last $indelpos next $nextindelpos\n";
while (@lines) {
my $line = pop(@lines);
$line =~ /^(.+?)\t(.+?)\t(.+?)\t(.+?)\t(.+?)\t/;
my $c = $1;
my $p = $2;
my $t = $3;
my $r = $4;
my $a = $5;
# print indels
if (length($r) - length($a) != 0) {
print $line;
} else {
# print other events which are more than prox away from indels
if (abs($indelpos - $p) >= $prox and abs($nextindelpos - $p) >= $prox) {
print $line;
}
}
}
$indelpos = $pos;
}
}
# flush lines end of file
while ($lines) {
print pop(@lines);
}
|