0
|
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 }
|