Mercurial > repos > ryanmorin > nextgen_variant_identification
comparison SNV/filter_snvmix.pl @ 0:74f5ea818cea
Uploaded
author | ryanmorin |
---|---|
date | Wed, 12 Oct 2011 19:50:38 -0400 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:74f5ea818cea |
---|---|
1 #!/usr/bin/perl | |
2 use strict; | |
3 use Getopt::Std; | |
4 use vars qw($opt_d $opt_p $opt_c $opt_C $opt_i $opt_S $opt_P $opt_I); | |
5 getopts("d:cC:i:p:SPI:"); | |
6 | |
7 my $usage = "$0 -p nonref_p-value_threshold [0.99] -c track_sequenced_bases -d minimum_nonref_depth [2] -C minimum_number_of_read_centered_calls [1] -i max_num_indels [0] -S need_dual_strand_coverage [0] -P use_unpaired_reads_too"; | |
8 | |
9 #chr:pos ref nref ref:ref_count,nref:nref_count,pAA,pAB,pBB,maxPr+ r- n+ n- nC nE nD | |
10 #chr:pos ref nref ref:ref_count,nref:nref_count,pAA,pAB,pBB,maxP ref_fwd ref_rev nref_fwd nref_rev nref_center nref_edges indel_pos ref_clean nref_clean ref_bad_pair nref_bad_pair ref_ins nref_ins ref_del nref_del ref_junc nref_junc nref_in_unique_pos | |
11 my $nonref_depth_thresh = $opt_d || 2; | |
12 my $min_centered = $opt_C || 1; | |
13 my $max_indels = $opt_i || 0; | |
14 my $dual_strand = $opt_S; | |
15 my $p_threshold = $opt_p || 0.99; | |
16 my $track_sequenced_bases = $opt_c; | |
17 my $pair_filtering = 1; | |
18 my $max_indel_relative = $opt_I; #filter on proportion of indels at SNV position relative to nonref calls. E.g. 0.5 would allow 1/2 the number of indel reads at that position vs indels. | |
19 if($opt_P){ | |
20 $pair_filtering = 0; | |
21 } | |
22 while(<STDIN>){ | |
23 chomp; | |
24 my ($chrpos,$ref,$nonref,$info,$ref_plus,$ref_minus,$nonref_plus,$nonref_minus,$centered,$ended,$indels,$ref_clean,$nonref_clean,$ref_bad_pair,$nonref_bad_pair,$ref_ins,$nonref_ins,$ref_del,$nonref_del,$ref_junc,$nref_junc,$nonref_in_unique_pos) = split; | |
25 my $ok; | |
26 my ($ref_info, $nref_info, $pAA, $pAB, $pBB, $call) = split(/,/, $info); | |
27 my $skip; | |
28 if($track_sequenced_bases){ | |
29 $skip = 1 if $call == 1; | |
30 } | |
31 else{ | |
32 next if $call == 1; | |
33 } | |
34 my $is_sequenced; | |
35 my ($foo,$nref_num) = split /:/, $nref_info; | |
36 if($pAA < ($pAB + $pBB)) { | |
37 if( ($pAB + $pBB) >= $p_threshold) { | |
38 $is_sequenced = $pAB+$pBB; | |
39 } | |
40 else{ | |
41 next; | |
42 } | |
43 } | |
44 else{ | |
45 if($track_sequenced_bases){ | |
46 if($pAA >= $p_threshold){ | |
47 $is_sequenced = $pAA; | |
48 } | |
49 # print STDERR "skip $skip\n"; | |
50 $skip = 1; | |
51 } | |
52 else{ | |
53 next; | |
54 } | |
55 | |
56 } | |
57 | |
58 if($track_sequenced_bases){ | |
59 if($is_sequenced){ | |
60 print STDERR "$chrpos\t$is_sequenced\n"; | |
61 } | |
62 next if $skip; | |
63 } | |
64 if($indels <= $max_indels && $centered >= $min_centered && $nref_num >= $nonref_depth_thresh){ | |
65 if($dual_strand){ | |
66 $ok=1 if $nonref_minus > 0 && $nonref_plus > 0; | |
67 } | |
68 else{ | |
69 $ok= 1; | |
70 } | |
71 } | |
72 else{ | |
73 #print "$_ FAIL\n$indels <= $max_indels && $centered >= $min_centered && $nref_num >= $nonref_depth_thresh\n"; | |
74 } | |
75 # if($nonref_bad_pair > ($nonref_plus+$nonref_minus)/2){ | |
76 if($nonref_bad_pair > 1 && $pair_filtering){ | |
77 #Rodrigo's bad pair filter fail | |
78 #print STDERR "skipping $chrpos due to $nonref_bad_pair > ($nonref_plus+$nonref_minus)/2\n"; | |
79 $ok = 0; | |
80 | |
81 } | |
82 | |
83 unless($chrpos =~ /^chr/){ | |
84 $chrpos = "chr$chrpos"; | |
85 } | |
86 if($ok){ | |
87 print "$chrpos\t$ref\t$nonref\t$info\n"; | |
88 } | |
89 } |