Mercurial > repos > mcharles > rapsosnp

comparison rapsodyn/mpileupfilterandstat.pl @ 7:3f7b0788a1c4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author mcharles
date Tue, 07 Oct 2014 10:34:34 -0400
parents 442a7c88b886
children
comparison
equal deleted inserted replaced
6:1776b8ddd87e 7:3f7b0788a1c4
1 #!/usr/bin/perl 1 #!/usr/bin/perl
2 #V1.0.0
2 use strict; 3 use strict;
3 use Getopt::Long; 4 use Getopt::Long;
4 5
5 # 6 #
6 # Filter a pileup file on forward/reverse presence and %read having the variant 7 # Filter a pileup file on forward/reverse presence and %read having the variant
26 my $MIN_DEPTH=0; 27 my $MIN_DEPTH=0;
27 my $MAX_DEPTH=500; 28 my $MAX_DEPTH=500;
28 my $VERBOSE=0; 29 my $VERBOSE=0;
29 my $ONLY_UNFILTERED_VARIANT="OFF"; 30 my $ONLY_UNFILTERED_VARIANT="OFF";
30 my $DO_STAT="NO"; 31 my $DO_STAT="NO";
32
33 my $nb_variant_checked=0;
34 my $nb_variant_selected=0;
31 35
32 36
33 my $STAT_MIN_DEPTH_MIN = 2; 37 my $STAT_MIN_DEPTH_MIN = 2;
34 my $STAT_MIN_DEPTH_MAX = 10; 38 my $STAT_MIN_DEPTH_MAX = 10;
35 my $STAT_MIN_DEPTH_STEP = 2; 39 my $STAT_MIN_DEPTH_STEP = 2;
81 85
82 86
83 #Extraction des variants 87 #Extraction des variants
84 my $nb_line=0; 88 my $nb_line=0;
85 while (my $line=<IF>){ 89 while (my $line=<IF>){
90 $nb_variant_checked++;
86 $nb_line++; 91 $nb_line++;
87 if (($nb_line % 1000000 == 0)&&($VERBOSE==1)){ 92 if (($nb_line % 1000000 == 0)&&($VERBOSE==1)){
88 print "$nb_line\n"; 93 print "$nb_line\n";
89 } 94 }
90 my $error_code=0; 95 my $error_code=0;
149 $after = $tbl_line[$i+1]; 154 $after = $tbl_line[$i+1];
150 } 155 }
151 my $error_code = check_error($tbl_line[$i],$before,$after,\%USR_PARAM); 156 my $error_code = check_error($tbl_line[$i],$before,$after,\%USR_PARAM);
152 if ($error_code == 0){ 157 if ($error_code == 0){
153 print $line{"line"}; 158 print $line{"line"};
159 $nb_variant_selected++;
154 } 160 }
155 else { 161 else {
156 push(@error,$error_code,"\t",$line{"line"}); 162 push(@error,$error_code,"\t",$line{"line"});
157 } 163 }
158 } 164 }
159 165
160 ### LOG 166 ### LOG
161 open(LF,">$logfile") or die ("Can't open $logfile\n"); 167 open(LF,">$logfile") or die ("Can't open $logfile\n");
162 168 print LF "\n####\t MPileup filtering \n";
163 if ($DO_STAT eq "YES"){ 169 print LF "Variant checked :\t$nb_variant_checked\n";
170 if ($DO_STAT eq "NO"){
171 print LF "Variant selected :\t$nb_variant_selected\n";
172 }
173 elsif ($DO_STAT eq "YES"){
164 for (my $idx_min_depth=$STAT_MIN_DEPTH_MIN;$idx_min_depth<=$STAT_MIN_DEPTH_MAX;$idx_min_depth = $idx_min_depth + $STAT_MIN_DEPTH_STEP ){ 174 for (my $idx_min_depth=$STAT_MIN_DEPTH_MIN;$idx_min_depth<=$STAT_MIN_DEPTH_MAX;$idx_min_depth = $idx_min_depth + $STAT_MIN_DEPTH_STEP ){
165 for (my $idx_max_depth=$STAT_MAX_DEPTH_MIN;$idx_max_depth<=$STAT_MAX_DEPTH_MAX;$idx_max_depth = $idx_max_depth + $STAT_MAX_DEPTH_STEP ){ 175 for (my $idx_max_depth=$STAT_MAX_DEPTH_MIN;$idx_max_depth<=$STAT_MAX_DEPTH_MAX;$idx_max_depth = $idx_max_depth + $STAT_MAX_DEPTH_STEP ){
166 for (my $idx_freq = $STAT_FREQ_MIN;$idx_freq<=$STAT_FREQ_MAX;$idx_freq= $idx_freq+$STAT_FREQ_STEP){ 176 for (my $idx_freq = $STAT_FREQ_MIN;$idx_freq<=$STAT_FREQ_MAX;$idx_freq= $idx_freq+$STAT_FREQ_STEP){
167 for (my $idx_dist=$STAT_DIST_MIN;$idx_dist<=$STAT_DIST_MAX;$idx_dist = $idx_dist + $STAT_DIST_STEP){ 177 for (my $idx_dist=$STAT_DIST_MIN;$idx_dist<=$STAT_DIST_MAX;$idx_dist = $idx_dist + $STAT_DIST_STEP){
168 for (my $idx_fr=0;$idx_fr<=1;$idx_fr++){ 178 for (my $idx_fr=0;$idx_fr<=1;$idx_fr++){
181 } 191 }
182 } 192 }
183 } 193 }
184 194
185 195
186 for (my $i=0;$i<=$#error;$i++){ 196 #for (my $i=0;$i<=$#error;$i++){
187 print LF $error[$i]; 197 # print LF $error[$i];
188 } 198 #}
189 close (LF); 199 close (LF);
190 200
191 201
192 202
193 203