#!/usr/bin/perl -w
use strict;

unless($ARGV[1]){
print <<EOF;
	perl define.sensitive.pl User_input Polymorphism_data(with Derived Allele Frequency)

EOF
exit;
}


my $user_annotation = $ARGV[0];
my $snp_data = $ARGV[1];
my %saw_snp;
my $daf_cut_off = 0.005;
my $r_exc = "a.R";
my $sensitive_out = "sensitive.user.bed";


&read_snp();
open(SEN,">$sensitive_out")||die;


my @samples = split /\n/, `cut -f 4 $user_annotation | sort|uniq`;
if (scalar @samples <=1){
	my $anno = `sort -k 1,1 -k 2,2g $user_annotation | mergeBed`;
	&cal($anno);
}else{
	foreach my $sample(@samples){
		my $anno = `awk '{FS="\t";OFS="\t"} \$4 == "$sample"' $user_annotation |sort -k 1,1 -k 2,2g | mergeBed`;
       	 	&cal($anno);
	}
}


sub read_snp{
	open(IN,$snp_data)||die;
	while(<IN>){
		my ($chr,$pos,$daf) = (split /\t/,$_)[0,1,3];
		$saw_snp{$chr}{$pos} = $daf;
	}
	close IN;
}

sub cal{
	my ($cate) = @_;
	my $snp_count = 0;
	my $snp_rare = 0;
	my @interval = split /\n/,$cate;
	foreach my $line (@interval){
		my ($chr,$start,$end) = (split /\t/,$line)[0..2];
		$chr =~ s/^chr//i;
		$chr =~ s/^/chr/;
		foreach my $pos($start .. $end-1){
			if (defined $saw_snp{$chr}{$pos}){
				$snp_count ++;
				if ($saw_snp{$chr}{$pos} < $daf_cut_off){
					$snp_rare ++;
				}
			}
		} 
	}
	close IN;

	if ($snp_count > 0 && $snp_rare >0){
		open(R,">$r_exc");
		print R "p.value = binom.test($snp_rare,$snp_count,p=15835777/26349801",',alternative="g")$p.value',"\n";
		print R "if (p.value <= 1e-5 && $snp_rare/$snp_count >= 0.6539634){\n";
		print R "print ('USEN')\n";
		print R "} else if (p.value <= 1e-3 && $snp_rare/$snp_count >= 0.6230845){\n";
		print R "print ('SEN')}\n";
		close R;
		my $ind = `Rscript $r_exc`;
		if ($ind =~ /USEN/){
			foreach my $line (@interval){
				print SEN join("\t",(split /\t/,$line)[0..2],"Ultra-Sensitive"),"\n";
			}
		}elsif($ind =~ /SEN/){
			foreach my $line (@interval){
                                print SEN join("\t",(split /\t/,$line)[0..2],"Sensitive"),"\n";
                        }
		}
	}
}


unlink($r_exc);
#if (-z $sensitive_out){
#	unlink($sensitive_out);
#}