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

Uploaded

24:e8e6b962c1f2

my $MAX_DEPTH=500;
my $VERBOSE=0;
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);
}

GetOptions (
"input_file=s" => \$inputfile,
"log_file=s" => \$logfile,
"min_depth=i" => \$MIN_DEPTH,

"min_frequency=f" => \$MIN_VARIANTFREQUENCY,
"min_distance=i" => \$MIN_DISTANCE,
"min_forward_and_reverse=i" => \$MIN_FORWARDREVERSE,
"variant_only=s" => \$ONLY_UNFILTERED_VARIANT,
"v=i" => \$VERBOSE,
"do_stat=s" => \$DO_STAT
) or die("Error in command line arguments\n");


open(IF, $inputfile)  or die("Can't open $inputfile\n");

my @tbl_line;
my %USR_PARAM;

### LOG
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;

					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";
	}
}


for (my $i=0;$i<=$#error;$i++){

25:39376c7204be

my $MAX_DEPTH=500;
my $VERBOSE=0;
my $ONLY_UNFILTERED_VARIANT="OFF";
my $DO_STAT="NO";


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,
"log_file=s" => \$logfile,
"min_depth=i" => \$MIN_DEPTH,

"min_frequency=f" => \$MIN_VARIANTFREQUENCY,
"min_distance=i" => \$MIN_DISTANCE,
"min_forward_and_reverse=i" => \$MIN_FORWARDREVERSE,
"variant_only=s" => \$ONLY_UNFILTERED_VARIANT,
"v=i" => \$VERBOSE,
"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;
my %USR_PARAM;

### LOG
open(LF,">$logfile") or die ("Can't open $logfile\n");

if ($DO_STAT eq "YES"){
	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 "\n";
		}
	}
}


for (my $i=0;$i<=$#error;$i++){