#!/usr/local/bin/perl 

use strict;
use Getopt::Long;

# Works best if you use -z 0.9 option when you do mapping. Assumes read length of 0.75. 

my $distance = 200;
my $threshold = 0;
my $bcf_file = '';
my $help = 0;
my $readlength = 100;

GetOptions ( 'i=s' => \$bcf_file, 
			 'd=s' =>  \$distance,
			 't=s' =>\$threshold,
			 'h' => \$help,
			 'r=s' =>\$readlength
				);
				
if ($help) { &help; }
#if ($bcf_file = 0) { &help; }

#help output
sub help {
	print "het_finder.pl identifies and lists possible mixed alleles/minority variants.\nOnly those of high mapping quality, and lacking strand bias will be included in an attempt to exclude those that may be due to sequencing or mapping error\n
	-i		input bcf file (_variant.bcf file will be less stringent but faster)
	-r		read length (default 100)
	-d 		variants within x base-pairs of one another will be filtered to remove likely mismapped regions (default 200)
	-t		Only include variants above threshold (between 0-1). Lower limit is usually ~0.04 anyway (default 0)\n";

	exit;
}


my $reference = '';
my $count = 0;

#get sanger ID from bcf_file
my $bcf = '';
if ($bcf_file =~ m/(\w*#\w*|\w*)\.bcf/){
#if ($bcf_file =~ m/\/(\w*)_variant.bcf/){
	 $bcf = "$1";
	 print "$bcf\n";
}

my $length = '';
my $to_print = '';
my $filename = '';
my @list_of_pos;
my @pos = ();
my @qual = ();
my @alt_all = ();
my @ref_all = ();
my @DP4 = ();
my $switch = 0;
my @PV4 = ();
my $dec = '';
my $ref = '';
my $alt = '';
my $cov = '';
my $maxcov = 500;
my $mincov = 10;
my $percent = '';
my @base2 = ();
my @base1 = ();
my $bam = '';
my $mpileup_out = '';
my $string = '';
my ($A_count, $T_count,  $G_count, $C_count, $Aa_count,$Tt_count, $Cc_count, $Gg_count, $g_count, $c_count, $t_count, $a_count) = '';
my ($lower, $higher) = '';
my $strandbias = 0;
my $mapQ = 0;
my $baseQ = 0;
my $endbias = 0;
my $upperthreshold = (1 - $threshold);

$to_print .= "Position\tDP4_in_bcf\tsupport_in_bam:A,T,C,Gs\tfrequency_of_ALT\n";

#Get expected coverage and set -d and -D for vcfutils - assumes readlength of 75

#if ($bcf_file =~ m/(\w*)\.bcf/){
if ($bcf_file =~ m/^(.*).bcf/){
	 $bam = "$1";
#	 print "$bam";
}
if ($bcf_file =~ m/^(.*)_variant.bcf/){
	 $bam = "$1";
	 print "$bam\n";
}

my $stats = `samtools idxstats $bam.bam`;
if ($stats =~ m/^.*\t(\d+)\t(\d+)\t/){
	$cov=(($2*$readlength)/$1);
	$mincov=($cov/2);
	$maxcov=($cov*2);
	}


#use bcftools & vcfutils to extract file considering coverage and strand bias etc (simon has min P value of 0.001 for strand bias)

	system("bcftools view  $bcf_file | ~jb31/scripts/vcfutils_josie.pl varFilter -d 4   >> $bcf.temp");
	$filename = ("$bcf.temp");
	open( MIDFILE, "<$filename" ) or die "$!";
	my	@bcf_to_pars = <MIDFILE>;
	close (MIDFILE);

## fix - to account for multiple references just take the first one (eg plasmids)
if ($bcf_to_pars[50] =~ m/^(\S+)\t(\d+).*\t(\d+)\t.\tDP/){
$reference = "$1";
print "$reference\n";
}
	
#pars out hets and relevant info
	
	for my $k (0..$#bcf_to_pars){
		chomp $bcf_to_pars[$k];
		
#get rid of pos with 100% agreement with ref or ALT
		next if ($bcf_to_pars[$k] =~ m/DP4=\d+,\d+,0,0;/);
		next if ($bcf_to_pars[$k] =~ m/DP4=0,0,\d+,\d+;/);
		
#pull out info of other pos

		if ($bcf_to_pars[$k] =~ m/^(\S+)\t(\d+).*\t(\d+)\t.\tDP.*AF1=.*;AC.*DP4=(\d+,\d+,\d+,\d+).*PV4=(.*,.*,.*,.*)\s/){
			$reference = "$1";
			push (@pos, $2);
			push (@qual, $3);
			push (@DP4, $4);
			push (@PV4, $5);
							if ($4 =~ m/(\d+),(\d+),(\d+),(\d+)/){
							$ref = ($1+$2);
							push (@ref_all, $ref);
							$alt = ($3+$4);
							push (@alt_all, $alt);
							}
		}
	}


	
##################################
#filter for min coverage of 4 with atleast 1 on each strand (for both ref and alt) also can't 
#have any het positions within 200 bp (have to have more than 2 alt bases though to allow for errors) 
#also check for strand bias - have it at this point so can keep the adjacent snps in the array
	
	for my $k (0..$#pos){
			$switch = 0;
			#print "$pos[$k]";
			
			$lower = ($pos[$k]-$distance);
			$higher = ($pos[$k]+$distance);
					foreach my $m (0..$#pos){
					if (($pos[$m]< $pos[$k]) && ($pos[$m]> $lower) 
					&& ($alt_all[$m] > 4)
					) {
						$switch = 1;
						}

					if (($pos[$m]> $pos[$k]) && ($pos[$m]< $higher) 
				&& ($alt_all[$m] > 4)
					){
						$switch = 1;
						}
					}
			if ($PV4[$k] =~ m/(.*),(.*),(.*),(.*)/){
			$strandbias = $1;
			$baseQ = $2;
			$mapQ = $3;
			$endbias = $4;
			}
						
			if ($DP4[$k] =~ m/(\d+),(\d+),(\d+),(\d+)/){
			}		

			if (($ref_all[$k]>=4) && ($alt_all[$k]>=4) 
				&& ($1>1) && ($2>1) && ($3>1) && ($4>1) 
				&& ($qual[$k]> 100) 
				&& ($strandbias > 0.05) && ($baseQ > 0.001) && ($mapQ > 0.001) && ($endbias > 0.001) 
				&& (($ref_all[$k] + $alt_all[$k]) < $maxcov) && (($ref_all[$k] + $alt_all[$k]) > $mincov) 
				&& ($switch == 0)	) {
																						
																																						
			
			print "$lower $pos[$k-1] $higher $pos[$k+1]\n";
			
			$percent = sprintf("%.3f",($alt_all[$k]/($ref_all[$k]+$alt_all[$k])));
			
#get mpileup output and count number of bases			
			
			$mpileup_out = `samtools mpileup -r '$reference:$pos[$k]-$pos[$k]' $bam.bam`;
			
				if ($mpileup_out =~ m/\t$pos[$k]\t\w\t\d+\t(.*)\t/){
				$string = $1;
				$A_count = () = $string =~ m/A/g;
				$a_count = () = $string =~ m/a/g;
				$T_count = () = $string =~ m/T/g;
				$t_count = () = $string =~ m/t/g;
				$C_count = () = $string =~ m/C/g;
				$c_count = () = $string =~ m/c/g;
				$G_count = () = $string =~ m/G/g;
				$g_count = () = $string =~ m/g/g;
				
				$Aa_count = ($A_count + $a_count);
				$Tt_count = ($T_count + $t_count);
				$Cc_count = ($C_count + $c_count);
				$Gg_count = ($G_count + $g_count);

# Check that for any allele at frequency greater than 1 (so not the allowence of 1 error) is on both strands

					if ($Aa_count > 1 and ($A_count < 2 or $a_count < 2)){
					next;
					}
					if ($Tt_count > 1 and ($T_count < 2 or $t_count < 2)){
					next;
					}
					if ($Cc_count > 1 and ($C_count < 2 or $c_count < 2)){
					next;
					}
					if ($Gg_count > 1 and ($G_count < 2 or $g_count < 2)){
					next;
					}
				
# Makesure there isn't more than 1 read covering bases that arn't ALT or REF or there arnt all bases
 					if (($Aa_count > 1 and $Tt_count > 1 and $Cc_count > 1) or ($Aa_count > 1 and $Tt_count > 1 and $Gg_count > 1)
 					or ($Aa_count > 1 and $Cc_count > 1 and $Gg_count > 1) or ($Tt_count > 1 and $Cc_count > 1 and $Gg_count > 1)
 					or ($Tt_count > 0 and $Cc_count > 0 and $Gg_count > 0 and $Aa_count > 0)){
 					next;
 					}
 					elsif ($percent >= $threshold && $percent <= $upperthreshold ) { 
					$to_print .= "$pos[$k]\t$DP4[$k]\t$Aa_count,$Tt_count,$Cc_count,$Gg_count\t$percent\t$bcf\n"; 
 					}
				}
			
			}
			
			

	}		
		
	open NEWFILE, ">$bcf\_het\_search.out" or die "Can't write to file [$!]\n";
	print NEWFILE "$to_print";
	close NEWFILE;
	
	
	system("rm $bcf.temp");
	




