#!/usr/bin/perl

use strict;


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


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

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

my (%fr, %pair, %wild);


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



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

    my $fname = join('', 'alleles/alleles.17strains.chr', $ch, '.txt');


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

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

    $line = <IN>;
    my $header_line = $line;
    chomp($header_line);
    print("$header_line\n");
    my @header = split("\t", $header_line);
    print("$header[19]\n");

    while ( defined($line = <IN>) ) {
	chomp($line);
	print("$line\n");

        @row = split("\t", $line);
	$alt = $row[3];

	if ( length($alt) > 1 ) {
	    print("alt = $alt\n");
	    exit;
	} else {
	    if ( $row[19] eq '.' ) {
		print("spret allele missing: $row[19] - skip this marker\n");
	    } else {

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

		my @recode;

		if ( $row[19] == 0 ) {

		    $spret = $ref;
		    $der = $alt;
		    @recode = @row[4..20];

		} elsif ( $row[19] == 1 ) {

		    $spret = $alt;
		    $der = $ref;
		    my @oldall = @row[4..20];

		    for ($i=0; $i<17; $i++) {
			if ( $oldall[$i] eq '.' ) {
			    push(@recode, '.');
			} elsif ( $oldall[$i] == 0 ) {
			    push(@recode, 1);
			} elsif ( $oldall[$i] == 1 ) {
			    push(@recode, 0);
			} else {
			    print("ERROR - unrecognised allele $oldall[$i]\n");
			    exit;
			}
		    }

		    print("recoded allele vector @recode\n");
		    exit;

		} else {
		    print("ERROR - unrecognised spret allele $row[19]\n");
		    exit;
		}

	    } # end of skip marker loop - spret allele missing
	} # end of error check loop - multiple ALT alleles

	exit;

    } # end of while loop

    exit;

    close(IN);

    exit;

} # end of chr loop




#    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");
