diff rapsodyn/mpileupfilterandstat.pl @ 25:39376c7204be draft

Uploaded

--- a/rapsodyn/mpileupfilterandstat.pl	Fri Sep 05 06:12:10 2014 -0400
+++ b/rapsodyn/mpileupfilterandstat.pl	Wed Sep 10 05:12:05 2014 -0400
@@ -29,19 +29,19 @@
 my $ONLY_UNFILTERED_VARIANT="OFF";
 my $DO_STAT="NO";
 
-if ($#ARGV<0){
-	print "\n";
-	print "perl 020_FilterPileupv6 -input_file <mpileup_file> [OPTION]\n";
-	print "-input_file \tinputfile in mpileup format\n";
-	print "-log_file \tlogfile containing discarded mpileup lines and the errorcode associated\n";
-	print "-min_depth \tminimum depth required [1]\n";
-	print "-max_depth \tmaximim depth (position with more coverage will be discarded) [100]\n";
-	print "-min_frequency \tminimum variant frequency (0->1) [1] (default 1 => 100% reads show the variant at this position)\n";
-	print "-min_distance \tminimum distance between variant [0]\n";
-	print "-min_forward_and_reverse \tminimum number of reads in forward and reverse covering the variant required [0]\n";
-	print "\n";
-	exit(0);
-}
+
+my $STAT_MIN_DEPTH_MIN = 2;
+my $STAT_MIN_DEPTH_MAX = 10;
+my $STAT_MIN_DEPTH_STEP = 2;
+my $STAT_MAX_DEPTH_MIN = 100;
+my $STAT_MAX_DEPTH_MAX = 200;
+my $STAT_MAX_DEPTH_STEP = 100;
+my $STAT_FREQ_MIN = 0.8;
+my $STAT_FREQ_MAX = 1;
+my $STAT_FREQ_STEP = 0.1;
+my $STAT_DIST_MIN = 0;
+my $STAT_DIST_MAX = 50;
+my $STAT_DIST_STEP = 50;
 
 GetOptions (
 "input_file=s" => \$inputfile,
@@ -53,10 +53,21 @@
 "min_forward_and_reverse=i" => \$MIN_FORWARDREVERSE,
 "variant_only=s" => \$ONLY_UNFILTERED_VARIANT,
 "v=i" => \$VERBOSE,
-"do_stat=s" => \$DO_STAT
+"do_stat=s" => \$DO_STAT,
+"stat_min_depth_min=i" => \$STAT_MIN_DEPTH_MIN,
+"stat_min_depth_max=i" => \$STAT_MIN_DEPTH_MAX,
+"stat_min_depth_step=i" => \$STAT_MIN_DEPTH_STEP,
+"stat_max_depth_min=i" => \$STAT_MAX_DEPTH_MIN,
+"stat_max_depth_max=i" => \$STAT_MAX_DEPTH_MAX,
+"stat_max_depth_step=i" => \$STAT_MAX_DEPTH_STEP,
+"stat_freq_min=f" => \$STAT_FREQ_MIN,
+"stat_freq_max=f" => \$STAT_FREQ_MAX,
+"stat_freq_step=f" => \$STAT_FREQ_STEP,
+"stat_dist_min=i" => \$STAT_DIST_MIN,
+"stat_dist_max=i" => \$STAT_DIST_MAX,
+"stat_dist_step=i" => \$STAT_DIST_STEP
 ) or die("Error in command line arguments\n");
 
-
 open(IF, $inputfile)  or die("Can't open $inputfile\n");
 
 my @tbl_line;
@@ -150,22 +161,24 @@
 open(LF,">$logfile") or die ("Can't open $logfile\n");
 
 if ($DO_STAT eq "YES"){
-	for (my $idx_min_depth=2;$idx_min_depth<=32;$idx_min_depth = $idx_min_depth*2){
-		for (my $idx_freq = 0.5;$idx_freq<=1;$idx_freq= $idx_freq+0.1){ 
-			for (my $idx_dist=0;$idx_dist<=50;$idx_dist = $idx_dist + 50){
-				for (my $idx_fr=0;$idx_fr<=1;$idx_fr++){
-					my %stat_param;
-					$stat_param{"min_depth"}=$idx_min_depth;
-					$stat_param{"max_depth"}=250;
-					$stat_param{"min_freq"}=$idx_freq;
-					$stat_param{"min_fr"}=$idx_fr;
-					$stat_param{"min_dist"}=$idx_dist;
+	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 ){
+		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 ){
+			for (my $idx_freq = $STAT_FREQ_MIN;$idx_freq<=$STAT_FREQ_MAX;$idx_freq= $idx_freq+$STAT_FREQ_STEP){ 
+				for (my $idx_dist=$STAT_DIST_MIN;$idx_dist<=$STAT_DIST_MAX;$idx_dist = $idx_dist + $STAT_DIST_STEP){
+					for (my $idx_fr=0;$idx_fr<=1;$idx_fr++){
+						my %stat_param;
+						$stat_param{"min_depth"}=$idx_min_depth;
+						$stat_param{"max_depth"}=$idx_max_depth;
+						$stat_param{"min_freq"}=$idx_freq;
+						$stat_param{"min_fr"}=$idx_fr;
+						$stat_param{"min_dist"}=$idx_dist;
 
-					print LF "#SNP = ",&test_check(\@tbl_line,\%stat_param),"\tdepth (min/max) = ",$stat_param{"min_depth"}," / ",$stat_param{"max_depth"},"\tmin_dist=",$stat_param{"min_dist"},"\tmin_freq=",$stat_param{"min_freq"},"\tmin_forwardreverse = ",$stat_param{"min_fr"},"\n";
-				}
-			}	
+						print LF "#SNP = ",&test_check(\@tbl_line,\%stat_param),"\tdepth (min/max) = ",$stat_param{"min_depth"}," / ",$stat_param{"max_depth"},"\tmin_dist=",$stat_param{"min_dist"},"\tmin_freq=",$stat_param{"min_freq"},"\tmin_forwardreverse = ",$stat_param{"min_fr"},"\n";
+					}
+				}	
+			}
+			print "\n";
 		}
-		print "\n";
 	}
 }