48 lines
946 B
Perl
48 lines
946 B
Perl
|
#!/usr/bin/env perl
|
||
|
|
|
||
|
|
use strict;
|
||
|
|
use warnings;
|
||
|
|
|
||
|
|
|
||
|
|
die "usage:perl $0 vcf1r vcf2r outvcf" unless @ARGV == 3;
|
||
|
|
my ($vcf1r, $vcf2r, $outvcf) = @ARGV;
|
||
|
|
|
||
|
|
# open R1, "${outputdir}/mutation/${tumor}.1r.snp.indel.vcf";
|
||
|
|
# open R2, "${outputdir}/mutation/${tumor}.2r.snp.indel.vcf";
|
||
|
|
# open OUT, ">${outputdir}/mutation/${tumor}.snp.indel.vcf";
|
||
|
|
|
||
|
|
open R1, "$vcf1r";
|
||
|
|
open R2, "$vcf2r";
|
||
|
|
open OUT, ">$outvcf";
|
||
|
|
|
||
|
|
my (%r1, %r2);
|
||
|
|
while (<R2>) {
|
||
|
|
next if /^#/;
|
||
|
|
chomp;
|
||
|
|
my @line = split;
|
||
|
|
my $key = join("_", @line[0, 1, 3, 4]);
|
||
|
|
$r2{$key} = $_;
|
||
|
|
}
|
||
|
|
|
||
|
|
while (<R1>) {
|
||
|
|
if (/^#/) {
|
||
|
|
print OUT;
|
||
|
|
next;
|
||
|
|
}
|
||
|
|
chomp;
|
||
|
|
my @line = split;
|
||
|
|
my $freq = (split(":", $line[-1]))[4];
|
||
|
|
my $dp = (split(":", $line[-1]))[1];
|
||
|
|
next if $dp < 50;
|
||
|
|
|
||
|
|
if ($freq < 0.01) {
|
||
|
|
my $key = join("_", @line[0, 1, 3, 4]);
|
||
|
|
if (exists $r2{$key}) {
|
||
|
|
print OUT "$_\n";
|
||
|
|
}
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
print OUT "$_\n";
|
||
|
|
}
|
||
|
|
}
|