annotate libexec/filterSomaticVariants.pl @ 7:3e8541ece3c7

trying to fix environment variable 0.0.1
author mini
date Fri, 26 Sep 2014 14:30:27 +0200
parents 87568e5a7d4f
children 0e8e6011082b
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
6
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
1 #!/usr/bin/env perl
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
2
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
3 =head1 LICENSE
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
4
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
5 Strelka Workflow Software
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
6 Copyright (c) 2009-2013 Illumina, Inc.
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
7
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
8 This software is provided under the terms and conditions of the
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
9 Illumina Open Source Software License 1.
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
10
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
11 You should have received a copy of the Illumina Open Source
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
12 Software License 1 along with this program. If not, see
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
13 <https://github.com/downloads/sequencing/licenses/>.
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
14
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
15 =head1 SYNOPSIS
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
16
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
17 filterSomaticVariants.pl [options] | --help
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
18
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
19 =head2 SUMMARY
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
20
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
21 Aggregate and filter the variant caller results by chromosome
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
22
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
23 =cut
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
24
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
25 use warnings FATAL => 'all';
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
26 use strict;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
27
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
28 use Carp;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
29 $SIG{__DIE__} = \&Carp::confess;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
30
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
31 use File::Spec;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
32 use Getopt::Long;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
33 use Pod::Usage;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
34
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
35 my $baseDir;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
36 my $libDir;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
37 #my $vcftDir;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
38 BEGIN {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
39 my $thisDir=(File::Spec->splitpath($0))[1];
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
40 $baseDir=File::Spec->catdir($thisDir,File::Spec->updir());
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
41 $libDir=File::Spec->catdir($baseDir,'lib');
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
42 #my $optDir=File::Spec->catdir($baseDir,'opt');
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
43 #$vcftDir=File::Spec->catdir($optDir,'vcftools','lib','perl5','site_perl');
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
44 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
45 use lib $libDir;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
46 use Utils;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
47 #use lib $vcftDir;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
48 #use Vcf;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
49
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
50
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
51 my $scriptName=(File::Spec->splitpath($0))[2];
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
52 my $argCount=scalar(@ARGV);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
53 my $cmdline = join(' ',$0,@ARGV);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
54
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
55
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
56 my ($chrom, $configFile);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
57 my $help;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
58
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
59 GetOptions( "chrom=s" => \$chrom,
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
60 "config=s" => \$configFile,
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
61 "help|h" => \$help) or pod2usage(2);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
62
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
63 pod2usage(2) if($help);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
64 pod2usage(2) unless(defined($chrom));
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
65 pod2usage(2) unless(defined($configFile));
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
66
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
67
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
68
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
69 #
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
70 # read config and validate values
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
71 #
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
72 checkFile($configFile,"configuration ini");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
73 my $config = parseConfigIni($configFile);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
74
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
75 my $chromSizeKey = "chrom_" . $chrom . "_size";
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
76 my $chromKnownSizeKey = "chrom_" . $chrom . "_knownSize";
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
77
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
78 for (("outDir", $chromSizeKey, $chromKnownSizeKey)) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
79 errorX("Undefined configuration option: '$_'") unless(defined($config->{derived}{$_}));
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
80 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
81 for (qw(binSize ssnvQuality_LowerBound sindelQuality_LowerBound
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
82 isSkipDepthFilters depthFilterMultiple
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
83 snvMaxFilteredBasecallFrac snvMaxSpanningDeletionFrac
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
84 indelMaxRefRepeat indelMaxWindowFilteredBasecallFrac
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
85 indelMaxIntHpolLength)) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
86 errorX("Undefined configuration option: '$_'") unless(defined($config->{user}{$_}));
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
87 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
88
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
89 my $outDir = $config->{derived}{outDir};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
90 my $chromDir = File::Spec->catdir($outDir,'chromosomes',$chrom);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
91 checkDir($outDir,"output");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
92 checkDir($chromDir,"output chromosome");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
93
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
94 my $userconfig = $config->{user};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
95
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
96 my $binList = getBinList($config->{derived}{$chromSizeKey},$userconfig->{binSize});
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
97 for my $binId (@$binList) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
98 my $dir = File::Spec->catdir($chromDir,'bins',$binId);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
99 checkDir($dir,"input bin");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
100 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
101
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
102 #
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
103 # parameters from user config file:
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
104 #
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
105
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
106 # minimum passed ssnv_nt Q-score:
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
107 my $minQSSNT=$userconfig->{ssnvQuality_LowerBound};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
108 # minimum passed sindel_nt Q-score:
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
109 my $minQSINT=$userconfig->{sindelQuality_LowerBound};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
110
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
111 #
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
112 # filtration parameters from user config file:
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
113 #
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
114
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
115 # skip depth filters for targeted resequencing:
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
116 my $isUseDepthFilter=(! $userconfig->{isSkipDepthFilters});
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
117 # multiple of the normal mean coverage to filter snvs and indels
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
118 my $depthFilterMultiple=$userconfig->{depthFilterMultiple};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
119 # max filtered basecall fraction for any sample
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
120 my $snvMaxFilteredBasecallFrac=$userconfig->{snvMaxFilteredBasecallFrac};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
121 # max snv spanning-deletion fraction for any sample
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
122 my $snvMaxSpanningDeletionFrac=$userconfig->{snvMaxSpanningDeletionFrac};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
123 # max indel reference repeat length
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
124 my $indelMaxRefRepeat=$userconfig->{indelMaxRefRepeat};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
125 # max indel window filter fraction
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
126 my $indelMaxWindowFilteredBasecallFrac=$userconfig->{indelMaxWindowFilteredBasecallFrac};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
127 # max indel interupted hompolymer length:
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
128 my $indelMaxIntHpolLength=$userconfig->{indelMaxIntHpolLength};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
129
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
130
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
131 # first we want the normal sample mean chromosome depth:
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
132 #
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
133 my $filterCoverage;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
134 if($isUseDepthFilter) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
135 my $totalCoverage = 0;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
136 for my $binId (@$binList) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
137 my $sFile = File::Spec->catfile($chromDir,'bins',$binId,'strelka.stats');
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
138 checkFile($sFile,"strelka bin stats");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
139 open(my $sFH, '<', $sFile)
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
140 || errorX("Can't open file: '$sFile' $!");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
141
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
142 my $is_found=0;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
143 while(<$sFH>) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
144 next unless(/^NORMAL_NO_REF_N_COVERAGE\s/);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
145 my $is_bad = 0;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
146 if(not /sample_size:\s*(\d+)\s+/) { $is_bad=1; }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
147 my $ss=$1;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
148 # leave the regex for mean fairly loose to pick up scientific notation, etc..
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
149 if(not /mean:\s*(\d[^\s]*|nan)\s+/) { $is_bad=1; }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
150 my $mean=$1;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
151 errorX("Unexpected format in file: '$sFile'") if($is_bad);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
152
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
153 my $coverage = ( $ss==0 ? 0 : int($ss*$mean) );
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
154
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
155 $totalCoverage += $coverage;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
156 $is_found=1;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
157 last;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
158 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
159 close($sFH);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
160 errorX("Unexpected format in file: '$sFile'") unless($is_found);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
161 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
162
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
163 my $chromKnownSize = $config->{derived}{$chromKnownSizeKey};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
164 my $normalMeanCoverage = ($totalCoverage/$chromKnownSize);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
165 $filterCoverage = $normalMeanCoverage*$depthFilterMultiple;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
166 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
167
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
168
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
169 # add filter description to vcf header unless it already exists
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
170 # return 1 if filter id already exists, client can decide if this is an error
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
171 #
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
172 sub add_vcf_filter($$$) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
173 my ($vcf,$id,$desc) = @_;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
174 return 1 if(scalar(@{$vcf->get_header_line(key=>'FILTER', ID=>$id)}));
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
175 $vcf->add_header_line({key=>'FILTER', ID=>$id,Description=>$desc});
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
176 return 0;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
177 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
178
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
179
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
180 sub check_vcf_for_sample($$$) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
181 my ($vcf,$sample,$file) = @_;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
182 my $is_found=0;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
183 for ($vcf->get_samples()) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
184 $is_found=1 if($_ eq $sample);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
185 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
186 errorX("Failed to find sample '$sample' in vcf file '$file'") unless($is_found);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
187 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
188
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
189
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
190
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
191 my $depthFiltId="DP";
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
192
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
193
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
194 # Runs all post-call vcf filters:
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
195 #
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
196 sub filterSnvFileList(\@$$$) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
197 my ($ifiles,$depthFilterVal,$acceptFileName,$isUseDepthFilter) = @_;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
198
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
199 my $baseFiltId="BCNoise";
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
200 my $spanFiltId="SpanDel";
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
201 my $qFiltId="QSS_ref";
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
202
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
203 open(my $aFH,'>',$acceptFileName)
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
204 or errorX("Failed to open file: '$acceptFileName'");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
205
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
206 my $is_first=1;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
207 for my $ifile (@$ifiles) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
208 open(my $iFH,'<',"$ifile")
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
209 or errorX("Failed to open file: '$ifile'");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
210 my $vcf = Vcf->new(fh=>$iFH);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
211
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
212 # run some simple header validation on each vcf:
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
213 $vcf->parse_header();
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
214 check_vcf_for_sample($vcf,'NORMAL',$ifile);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
215 check_vcf_for_sample($vcf,'TUMOR',$ifile);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
216
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
217 if($is_first) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
218 # TODO: update vcf meta-data for chromosome-level filtered files
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
219 #
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
220 $vcf->remove_header_line(key=>"cmdline");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
221 $vcf->add_header_line({key=>"cmdline",value=>$cmdline});
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
222 if($isUseDepthFilter) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
223 $vcf->add_header_line({key=>"maxDepth_$chrom",value=>$depthFilterVal});
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
224 add_vcf_filter($vcf,$depthFiltId,"Greater than ${depthFilterMultiple}x chromosomal mean depth in Normal sample");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
225 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
226 add_vcf_filter($vcf,$baseFiltId,"Fraction of basecalls filtered at this site in either sample is at or above $snvMaxFilteredBasecallFrac");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
227 add_vcf_filter($vcf,$spanFiltId,"Fraction of reads crossing site with spanning deletions in either sample exceeeds $snvMaxSpanningDeletionFrac");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
228 add_vcf_filter($vcf,$qFiltId,"Normal sample is not homozygous ref or ssnv Q-score < $minQSSNT, ie calls with NT!=ref or QSS_NT < $minQSSNT");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
229 print $aFH $vcf->format_header();
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
230 $is_first=0;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
231 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
232
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
233 while(my $x=$vcf->next_data_hash()) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
234 my $norm=$x->{gtypes}->{NORMAL};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
235 my $tumr=$x->{gtypes}->{TUMOR};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
236
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
237 my %filters;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
238
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
239 # normal depth filter:
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
240 my $normalDP=$norm->{DP};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
241 if($isUseDepthFilter) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
242 $filters{$depthFiltId} = ($normalDP > $depthFilterVal);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
243 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
244
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
245 # filtered basecall fraction:
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
246 my $normal_filt=($normalDP>0 ? $norm->{FDP}/$normalDP : 0);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
247
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
248 my $tumorDP=$tumr->{DP};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
249 my $tumor_filt=($tumorDP>0 ? $tumr->{FDP}/$tumorDP : 0);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
250
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
251 $filters{$baseFiltId}=(($normal_filt >= $snvMaxFilteredBasecallFrac) or
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
252 ($tumor_filt >= $snvMaxFilteredBasecallFrac));
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
253
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
254 # spanning deletion fraction:
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
255 my $normalSDP=$norm->{SDP};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
256 my $normalSpanTot=($normalDP+$normalSDP);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
257 my $normalSpanDelFrac=($normalSpanTot>0 ? ($normalSDP/$normalSpanTot) : 0);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
258
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
259 my $tumorSDP=$tumr->{SDP};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
260 my $tumorSpanTot=($tumorDP+$tumorSDP);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
261 my $tumorSpanDelFrac=($tumorSpanTot>0 ? ($tumorSDP/$tumorSpanTot) : 0);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
262
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
263 $filters{$spanFiltId}=(($normalSpanDelFrac > $snvMaxSpanningDeletionFrac) or
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
264 ($tumorSpanDelFrac > $snvMaxSpanningDeletionFrac));
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
265
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
266 # Q-val filter:
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
267 $filters{$qFiltId}=(($x->{INFO}->{NT} ne "ref") or
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
268 ($x->{INFO}->{QSS_NT} < $minQSSNT));
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
269
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
270 $x->{FILTER} = $vcf->add_filter($x->{FILTER},%filters);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
271
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
272 print $aFH $vcf->format_line($x);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
273 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
274
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
275 $vcf->close();
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
276 close($iFH);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
277 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
278
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
279 close($aFH);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
280 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
281
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
282
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
283
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
284 sub updateA2(\@$) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
285 my ($a2,$FH) = @_;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
286 my $line=<$FH>;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
287 unless(defined($line)) { errorX("Unexpected end of somatic indel window file"); }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
288 chomp $line;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
289 @$a2 = split("\t",$line);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
290 unless(scalar(@$a2)) { errorX("Unexpected format in somatic indel window file"); }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
291 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
292
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
293
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
294
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
295 sub filterIndelFileList(\@$$$) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
296 my ($ifiles,$depthFilterVal,$acceptFileName,$isUseDepthFilter) = @_;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
297
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
298 my $repeatFiltId="Repeat";
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
299 my $iHpolFiltId="iHpol";
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
300 my $baseFiltId="BCNoise";
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
301 my $qFiltId="QSI_ref";
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
302
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
303 open(my $aFH,'>',$acceptFileName)
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
304 or errorX("Failed to open file: '$acceptFileName'");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
305
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
306 my $is_first=1;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
307 for my $ifile (@$ifiles) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
308 open(my $iFH,'<',"$ifile")
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
309 or errorX("Failed to open somatic indel file: '$ifile'");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
310
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
311 my $iwfile = $ifile . ".window";
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
312 open(my $iwFH,'<',"$iwfile")
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
313 or errorX("Failed to open somatic indel window file: '$iwfile'");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
314
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
315 my @a2; # hold window file data for one line in case we overstep...
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
316
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
317 my $vcf = Vcf->new(fh=>$iFH);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
318
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
319 # run some simple header validation on each vcf:
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
320 $vcf->parse_header();
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
321 check_vcf_for_sample($vcf,'NORMAL',$ifile);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
322 check_vcf_for_sample($vcf,'TUMOR',$ifile);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
323
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
324 if($is_first) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
325 # TODO: update all vcf meta-data for chromosome-level filtered files
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
326 #
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
327 $vcf->remove_header_line(key=>"cmdline");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
328 $vcf->add_header_line({key=>"cmdline",value=>$cmdline});
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
329 if($isUseDepthFilter) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
330 $vcf->add_header_line({key=>"maxDepth_$chrom",value=>$depthFilterVal});
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
331 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
332 $vcf->add_header_line({key=>'FORMAT', ID=>'DP50',Number=>1,Type=>'Float',Description=>'Average tier1 read depth within 50 bases'});
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
333 $vcf->add_header_line({key=>'FORMAT', ID=>'FDP50',Number=>1,Type=>'Float',Description=>'Average tier1 number of basecalls filtered from original read depth within 50 bases'});
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
334 $vcf->add_header_line({key=>'FORMAT', ID=>'SUBDP50',Number=>1,Type=>'Float',Description=>'Average number of reads below tier1 mapping quality threshold aligned across sites within 50 bases'});
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
335 if($isUseDepthFilter) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
336 add_vcf_filter($vcf,$depthFiltId,"Greater than ${depthFilterMultiple}x chromosomal mean depth in Normal sample");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
337 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
338 add_vcf_filter($vcf,$repeatFiltId,"Sequence repeat of more than ${indelMaxRefRepeat}x in the reference sequence");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
339 add_vcf_filter($vcf,$iHpolFiltId,"Indel overlaps an interupted homopolymer longer than ${indelMaxIntHpolLength}x in the reference sequence");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
340 add_vcf_filter($vcf,$baseFiltId,"Average fraction of filtered basecalls within 50 bases of the indel exceeds $indelMaxWindowFilteredBasecallFrac");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
341 add_vcf_filter($vcf,$qFiltId,"Normal sample is not homozygous ref or sindel Q-score < $minQSINT, ie calls with NT!=ref or QSI_NT < $minQSINT");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
342 print $aFH $vcf->format_header();
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
343 $is_first=0;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
344 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
345
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
346 while(my $x=$vcf->next_data_hash()) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
347 $vcf->add_format_field($x,'DP50');
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
348 $vcf->add_format_field($x,'FDP50');
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
349 $vcf->add_format_field($x,'SUBDP50');
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
350
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
351 my $norm=$x->{gtypes}->{NORMAL};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
352 my $tumr=$x->{gtypes}->{TUMOR};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
353
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
354 my $chrom=$x->{CHROM};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
355 my $pos=int($x->{POS});
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
356
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
357 # get matching line from window file:
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
358 while((scalar(@a2)<2) or
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
359 (($a2[0] le $chrom) and (int($a2[1]) < $pos))) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
360 updateA2(@a2,$iwFH);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
361 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
362 unless(scalar(@a2) and ($a2[0] eq $chrom) and (int($a2[1]) == $pos))
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
363 { errorX("Can't find somatic indel window position.\nIndel line: " . $vcf->format_line($x) ); }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
364
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
365 # add window data to vcf record:
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
366 #
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
367 $norm->{DP50} = $a2[2]+$a2[3];
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
368 $norm->{FDP50} = $a2[3];
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
369 $norm->{SUBDP50} = $a2[4];
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
370 $tumr->{DP50} = $a2[5]+$a2[6];
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
371 $tumr->{FDP50} = $a2[6];
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
372 $tumr->{SUBDP50} = $a2[7];
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
373
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
374 my %filters;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
375
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
376 # normal depth filter:
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
377 my $normalDP=$norm->{DP};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
378 if($isUseDepthFilter) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
379 $filters{$depthFiltId}=($normalDP > $depthFilterVal);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
380 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
381
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
382 # ref repeat
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
383 my $refRep=$x->{INFO}->{RC};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
384 $filters{$repeatFiltId}=(defined($refRep) and
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
385 ($refRep > $indelMaxRefRepeat));
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
386
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
387 # ihpol
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
388 my $iHpol=$x->{INFO}->{IHP};
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
389 $filters{$iHpolFiltId}=(defined($iHpol) and
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
390 ($iHpol > $indelMaxIntHpolLength));
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
391
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
392 # base filt:
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
393 my $normWinFrac=( $norm->{DP50} ? $norm->{FDP50}/$norm->{DP50} : 0 );
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
394 my $tumrWinFrac=( $tumr->{DP50} ? $tumr->{FDP50}/$tumr->{DP50} : 0 );
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
395 $filters{$baseFiltId}=( ($normWinFrac >= $indelMaxWindowFilteredBasecallFrac) or
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
396 ($tumrWinFrac >= $indelMaxWindowFilteredBasecallFrac) );
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
397
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
398 # Q-val filter:
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
399 $filters{$qFiltId}=(($x->{INFO}->{NT} ne "ref") or
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
400 ($x->{INFO}->{QSI_NT} < $minQSINT));
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
401
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
402 $x->{FILTER} = $vcf->add_filter($x->{FILTER},%filters);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
403
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
404 print $aFH $vcf->format_line($x);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
405 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
406
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
407 $vcf->close();
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
408 close($iFH);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
409 close($iwFH);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
410 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
411
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
412 close($aFH);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
413 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
414
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
415
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
416
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
417 my @ssnvFiles;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
418 my @sindelFiles;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
419 for my $binId (@$binList) {
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
420 my $ssnvFile = File::Spec->catfile($chromDir,'bins',$binId,'somatic.snvs.unfiltered.vcf');
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
421 my $sindelFile = File::Spec->catfile($chromDir,'bins',$binId,'somatic.indels.unfiltered.vcf');
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
422 checkFile($ssnvFile,"bin snv file");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
423 checkFile($sindelFile,"bin indel file");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
424 push @ssnvFiles,$ssnvFile;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
425 push @sindelFiles,$sindelFile;
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
426 }
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
427
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
428 my $ssnvOutFile = File::Spec->catfile($chromDir,"somatic.snvs.vcf");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
429 filterSnvFileList(@ssnvFiles,$filterCoverage,$ssnvOutFile,$isUseDepthFilter);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
430
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
431 my $sindelOutFile = File::Spec->catfile($chromDir,"somatic.indels.vcf");
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
432 filterIndelFileList(@sindelFiles,$filterCoverage,$sindelOutFile,$isUseDepthFilter);
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
433
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
434
87568e5a7d4f Testing strelka version 0.0.1
mini
parents:
diff changeset
435 1;