#!/usr/bin/perl

use strict;


my @chrs = (1..19, 'X', 'Y');
#print("@chrs\n");

my $Ename = 'SNPs.17str.REFmismatches.txt';


my ($line, $ch, $i, $chr, $bppos);
my ($alt_freq, $ref, $geno_string);
my ($nskip, $nkeep, $nmatch, $ncomp, $nerr, $nwild, $nboth);

my (@row, @freq, @sev);

my (%fr, %pair);


# calculate freq and fr

for ($i=0; $i<=40; $i++) {
    $freq[$i] = $i/40;
}
#print("freq vector @freq\n");

@fr{ @freq } = ();
#foreach my $k ( keys %fr ) {
#    print("key $k\n");
#}

$pair{ A } = 'T';
$pair{ C } = 'G';
$pair{ G } = 'C';
$pair{ T } = 'A';



open(OUT, '>', $Ename) or die "Cannot open file $Ename for writing, reason $!\n";


foreach $ch ( @chrs ) {
    print("\nchromosome $ch\n\n");

    my $Wname = join('', 'stats/stats.known.chr', $ch, '.txt');
    my $Sname = join('', 'alleles/alleles.17strains.chr', $ch, '.txt');


    print("read in wild known data\n");

    my %wild;

    open(IN, '<', $Wname) or die "Cannot open file $Wname, reason $!\n";

    $line = <IN>;
    my $header_line = $line;
    chomp($header_line);
#    print("$header_line\n");

    $nskip = 0;
    while ( defined($line = <IN>) ) {

	chomp($line);
	@row = split(/\t/, $line);

	$alt_freq = $row[9];
#	print("alt freq = $alt_freq   ");
	if ( exists( $fr{ $alt_freq } ) ) {
#	    print("keep SNP\n");

	    $bppos = $row[1];
	    $wild{ $bppos }{'REF'} = $row[2];
	    $wild{ $bppos }{'ALT'} = $row[3];
	    $wild{ $bppos }{'ALT freq'} = $alt_freq;

	} else {
#	    print("missing data - awkward freq - drop SNP\n");
	    $nskip = ( $nskip + 1 );
	}
    } # end of while loop

    print("dropped $nskip SNPs\n\n");

    close(IN);


    print("read in 17 strains data\n");

    open(IN, '<', $Sname) or die "Cannot open file $Sname, reason $!\n";

    $line = <IN>;
    my $header_line = $line;
    chomp($header_line);
    my @header = split("\t", $header_line);
    my $error_header_line = join("\t", @header[0..1], 'WILDREF', @header[2..$#header]);
    if ( $ch == 1 ) {
	print(OUT "$error_header_line\n");
    }

    $nkeep = 0;
    $nmatch = 0;
    $ncomp = 0;
    $nerr = 0;
    while ( defined($line = <IN>) ) {

	chomp($line);
        @row = split("\t", $line);

	$bppos = $row[1];
	$ref = $row[2];

        if ( exists( $wild{ $bppos } ) ) {
	    # SNP exists in wild data
	    # check if REF matches

	    $nkeep = ( $nkeep + 1 );

	    if ( $ref eq ${$wild{ $bppos }}{'REF'} ) {

		$nmatch = ( $nmatch + 1 );

		$geno_string = join("\t", @row[4..$#row]);
		$wild{ $bppos }{'17str'} = $geno_string;

	    } elsif ( $pair{ $ref } eq ${$wild{ $bppos }}{'REF'} ) {
#		print("other strand of $ref is $pair{ $ref } which matches ${$wild{ $bppos }}{'REF'}\n");
		$ncomp = ( $ncomp + 1 );
		my $error_line = join("\t", @row[0..1], ${$wild{ $bppos }}{'REF'}, @row[2..$#row]);
		print(OUT "$error_line\n");

	    } else {
		$nerr = ( $nerr + 1 );
		my $error_line = join("\t", @row[0..1], ${$wild{ $bppos }}{'REF'}, @row[2..$#row]);
		print(OUT "$error_line\n");
#		print("REF doesn t match  ");
#		print("bp $bppos : wild ref = ${$wild{ $bppos }}{'REF'}  17 str ref = $ref\n");
	    }
	}

    } # end of while loop

    print("kept $nkeep SNPs\nof which $nmatch match REF\n$ncomp comp strand SNPs\nand $nerr REF errors\n\n");

    close(IN);


#    foreach my $k ( keys %fr ) {      
#	$fr{ $k }{'17 counts'} = [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ];
#	$fr{ $k }{'17 missing'} = [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ];
#	$fr{ $k }{'17 alts'} = [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ];
#    }


#    $nwild = 0;
#    foreach my $loc ( keys %wild ) {
#	$nwild = ( $nwild + 1 );
#	if ( exists( $wild{ $loc }{'17str'} )) {
#	    $nboth = ( $nboth + 1 );
#	    @sev = split("\t", ${$wild{ $loc }}{'17str'} );
#	    for ( $i=0; $i<=$#sev; $i++ ) {
#		if ( $sev[$i] eq '.' ) {
#		    # missing data
#		    $fr{ ${$wild{$loc}}{'ALT freq'} }{'17 missing'}[$i] = ( $fr{ ${$wild{$loc}}{'ALT freq'} }{'17 missing'}[$i] + 1 );
#		} elsif ( $sev[$i] == 0 || $sev[$i] == 1 ) {
#		    $fr{ ${$wild{$loc}}{'ALT freq'} }{'17 counts'}[$i] = ( $fr{ ${$wild{$loc}}{'ALT freq'} }{'17 counts'}[$i] + 1 );
#		    $fr{ ${$wild{$loc}}{'ALT freq'} }{'17 alts'}[$i] = ( $fr{ ${$wild{$loc}}{'ALT freq'} }{'17 alts'}[$i] + $sev[$i] );
#		} else {
#		    print("ERROR - $sev[$i] not recognised - should be ., 0 or 1\n");
#		}
#	    }
#	}
#    }

#    foreach my $k ( keys %fr ) {
#	for ( $i=0; $i<17; $i++) {
#	    print("i = $i   count = $fr{ $k }{'17 counts'}[$i]\n");
#	    $fr{ $k }{'17 alts'}[$i] = $fr{ $k }{'17 alts'}[$i]/$fr{ $k }{'17 counts'}[$i];
#	}
#    }

} # end of chr loop

close(OUT);



#    open(OUTK, '>', 'statistics.known.chr'.$ch.'.txt') or die "Cannot open file for writing, reason $!\n";                            
#    print(OUTK  "CHR\tBP\tREF\tALT\tMAF\tMAF_allele\tProp_missing_genotypes\tProp_missing_alleles\tAllfreq_ref\tAllfreq_alt\tGenofreq_ref\tGenofreq_het\tGenofreq_alt\tAllcount_ref\tAllcount_alt\tGenocount_ref\tGenocount_het\tGenocount_alt\n");

#       print(OUTK "$chr\t$bppos\t${$data{$bppos}}{'REF'}\t${$data{$bppos}}{'ALT'}\t${$data{$bppos}}{'MAF'}\t${$data{$bppos}}{'MAF allele'}\t${$data{$bppos}}{'missing'}[1]\t${$data{$bppos}}{'missing'}[0]\t${$data{$bppos}}{'allele freqs'}[0]\t${$data{$bppos}}{'allele freqs'}[1]\t${$data{$bppos}}{'genotype freqs'}[0]\t${$data{$bppos}}{'genotype freqs'}[1]\t${$data{$bppos}}{'genotype freqs'}[2]\t${$data{$bppos}}{'allele counts'}[0]\t${$data{$bppos}}{'allele counts'}[1]\t${$data{$bppos}}{'genotype counts'}[0]\t${$data{$bppos}}{'genotype counts'}[1]\t${$data{$bppos}}{'genotype counts'}[2]\n");
