-
Notifications
You must be signed in to change notification settings - Fork 1
/
SNP_GATK_Varscan_overlap.pl
executable file
·49 lines (43 loc) · 1.06 KB
/
SNP_GATK_Varscan_overlap.pl
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
#!/usr/bin/perl
use strict;
use warnings;
die "perl $0 <sample> <varscan> <gatk> <statfile>\n" if(@ARGV != 4);
my ($sample,$varscan,$gatk,$file)=@ARGV;
my %nsites;
my ($varsnp,$gatksnp,$gatkuniq,$varuniq,$overlap)=(0,0,0,0,0);
open IN1, ($varscan =~ /\.gz/? "gzip -dc $varscan |" : $varscan) || die $!;
<IN1>;
while(<IN1>){
chomp;
my @F = split /\t/;
my ($chr,$pos,$ref,$var)=@F[0,1,2,3];
my $id="$chr\t$pos\t$ref\t$var";
my $type = $F[12];
if($type eq "Germline"){
$nsites{$id}=1;
$varsnp++;
}
}
close IN1;
open IN2, ($gatk =~ /\.gz/? "gzip -dc $gatk |" : $gatk) || die $!;
chomp(my $head=<IN2>);
print "$head\n";
while(<IN2>){
chomp;
my @F = split /\t/;
my ($chr,$pos,$ref,$var)=@F[0,1,3,4];
my $id="$chr\t$pos\t$ref\t$var";
if(exists $nsites{$id}){
print "$_\n";
$overlap++;
}
else{
$gatkuniq++;
}
}
close IN2;
$gatksnp = $gatkuniq + $overlap;
$varuniq = $varsnp - $overlap;
open OUT,">>$file" || die $!;
print OUT "Samples\tGATK\tVarscan\tGATK_uniq\tVarscan_uniq\tOverlap\n";
print OUT "$sample\t$gatksnp\t$varsnp\t$gatkuniq\t$varuniq\t$overlap\n";