#!/usr/bin/perl

use strict;


my @chrs = (1..19);
#print("@chrs\n");
my @indices = (0..3, 6..25);


my ($fname, $line, $start, $ch);

my ($chr, $bppos, $mouse, $checksum, $den);
my ($count_ref, $count_alt, $count_mis, $ind_ref, $ind_alt, $ind_mis);
my ($count_hom_ref, $count_hom_alt, $count_het, $count_geno_mis, $geno_ref, $geno_het, $geno_alt, $geno_mis);
my ($freq_ref, $freq_alt, $freq_hom_ref, $freq_hom_alt, $freq_het, $prop_miss_all, $prop_miss_geno);
my (@row, @new_row, @mice, @mouse_info, @geno_vec);


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

    my $Kname = join('', 'snps.known.multALTs.chr', $ch, '.vcfstyle.txt');
    my $Nname = join('', 'snps.new.multALTs.chr', $ch, '.vcfstyle.txt');

    my $outKfile = 'stats/stats.known.chr'.$ch.'.txt';
    my $outNfile = 'stats/stats.new.chr'.$ch.'.txt';
    my $outAfile = 'stats/stats.chr'.$ch.'.txt';

    my %data;


    open(OUK, '>', $outKfile) or die "Cannot open file $outKfile for writing, reason $!\n";
    print(OUK  "CHR\tBP\tREF\tALT\tTRD\tProp_missing_genotypes\tProp_missing_alleles\tAllcount_ref\tGenocount_ref\n");

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

    $line = <IN>;
    $start = substr( $line, 0, 2);

    while ( $start eq '##' ) {
	$line = <IN>;
	$start = substr( $line, 0, 2);
    } # end of headers while loop

#    my $header_line = substr( $line, 1);
#    chomp($header_line);
#    my @new_header = @header[@indices];

    my $Nmult_knw = 0;

    while ( defined($line = <IN>) ) {

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

	$Nmult_knw++;

	@new_row = @row[@indices];

	$chr = $new_row[0];
	$bppos = $new_row[1];
	$data{ $bppos }{'new'} = 'N';
	$data{ $bppos }{'REF'} = $new_row[2];
	$data{ $bppos }{'ALT'} = $new_row[3];
	$data{ $bppos }{'TRD'} = $row[4];

	@mice = @new_row[ 4..$#new_row ];

	$count_ref = 0;
	$count_mis = 0;
	$count_hom_ref = 0;
	$count_geno_mis = 0;

	foreach $mouse ( @mice ) {
	    @mouse_info = split(':', $mouse);
	    @geno_vec = split("\/", $mouse_info[0]);
	    #	    print("genotype = @geno_vec\n");

	    $ind_ref = 0;
	    $ind_mis = 0;

	    if ( $geno_vec[0] eq '.' ) {
		$ind_mis++;
		$count_mis++;
	    } elsif ( $geno_vec[0] == 0 ) {
		$ind_ref = ( $ind_ref + 1 );
		$count_ref = ( $count_ref + 1 );
	    }

	    if ( $geno_vec[1] eq '.' ) {
		$ind_mis++;
		$count_mis++;
	    } elsif ( $geno_vec[1] == 0 ) {
		$ind_ref = ( $ind_ref + 1 );
		$count_ref = ( $count_ref + 1 );
	    }


	    $geno_mis = 0;
	    $geno_ref = 0;

	    if ( $ind_mis > 0 ) {
		$geno_mis++;
		$count_geno_mis++;
	    } elsif ( $ind_ref == 2 ) {
		$geno_ref++;
		$count_hom_ref++;
	    }

	} # end of loop over individual mice

	$data{ $bppos }{ 'allele counts' } = [ $count_ref, $count_mis ];
	$data{ $bppos }{ 'genotype counts' } = [ $count_hom_ref, $count_geno_mis ];

	$prop_miss_all = $count_mis/40;
	$prop_miss_geno = $count_geno_mis/20;
	$data{ $bppos }{ 'missing' } = [ $prop_miss_all, $prop_miss_geno ];

	print(OUK "$chr\t$bppos\t${$data{$bppos}}{'REF'}\t${$data{$bppos}}{'ALT'}\t${$data{$bppos}}{'TRD'}\t${$data{$bppos}}{'missing'}[1]\t${$data{$bppos}}{'missing'}[0]\t${$data{$bppos}}{'allele counts'}[0]\t${$data{$bppos}}{'genotype counts'}[0]\n");

    } # end of while loop
    close(IN);
    close(OUK);

    print("known stats file written  mult ALTs $Nmult_knw\n");


    open(OUN, '>', $outNfile) or die "Cannot open file $outNfile for writing, reason $!\n";
    print(OUN  "CHR\tBP\tREF\tALT\tTRD\tProp_missing_genotypes\tProp_missing_alleles\tAllcount_ref\tGenocount_ref\n");

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

    $line = <IN>;
    $start = substr( $line, 0, 2);

    while ( $start eq '##' ) {
      $line = <IN>;
      $start = substr( $line, 0, 2);
    }


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

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

      @new_row = @row[@indices];

      $chr = $new_row[0];
      $bppos = $new_row[1];
      $data{ $bppos }{'new'} = 'Y';
      $data{ $bppos }{'REF'} = $new_row[2];
      $data{ $bppos }{'ALT'} = $new_row[3];
      $data{ $bppos }{'TRD'} = $row[4];

      @mice = @new_row[ 4..$#new_row ];

      $count_ref = 0;
      $count_mis = 0;
      $count_hom_ref = 0;
      $count_geno_mis = 0;

      foreach $mouse ( @mice ) {

	  @mouse_info = split(':', $mouse);
	  @geno_vec = split("\/", $mouse_info[0]);
	  #	    print("genotype = @geno_vec\n");

	  $ind_ref = 0;
	  $ind_mis = 0;

	  if ( $geno_vec[0] eq '.' ) {
	      $ind_mis++;
	      $count_mis++;
	  } elsif ( $geno_vec[0] == 0 ) {
	      $ind_ref++;
	      $count_ref++;
	  }

	  if ( $geno_vec[1] eq '.' ) {
              $ind_mis++;
              $count_mis++;
          } elsif ( $geno_vec[1] == 0 ) {
              $ind_ref++;
              $count_ref++;
          }

	  $geno_mis = 0;
	  $geno_ref = 0;
	  if ( $ind_mis > 0 ) {
	      $geno_mis++;
	      $count_geno_mis++;
	  } elsif ( $ind_ref == 2 ) {
	      $geno_ref++;
	      $count_hom_ref++;
	  }
      } # end of loop over individual mice


      $data{ $bppos }{ 'allele counts' } = [ $count_ref, $count_mis ];
      $data{ $bppos }{ 'genotype counts' } = [ $count_hom_ref, $count_geno_mis ];

      $prop_miss_all = $count_mis/40;
      $prop_miss_geno = $count_geno_mis/20;
      $data{ $bppos }{ 'missing' } = [ $prop_miss_all, $prop_miss_geno ];

      print(OUN "$chr\t$bppos\t${$data{$bppos}}{'REF'}\t${$data{$bppos}}{'ALT'}\t${$data{$bppos}}{'TRD'}\t${$data{$bppos}}{'missing'}[1]\t${$data{$bppos}}{'missing'}[0]\t${$data{$bppos}}{'allele counts'}[0]\t${$data{$bppos}}{'genotype counts'}[0]\n");

    } # end of while loop
    close(IN);
    close(OUN);

    print("new stats file written mult ALTs $Nmult_new\n");


    my @srted_keys = sort { $a <=> $b } keys %data;

    open(OUA, '>', $outAfile) or die "Cannot open file $outAfile for writing, reason $!\n";
    print(OUA "CHR\tBP\tREF\tALT\tNEW\tTRD\tProp_missing_genotypes\tProp_missing_alleles\tAllcount_ref\tGenocount_ref\n");

    foreach my $loc ( @srted_keys ) {
	print(OUA "$chr\t$loc\t${$data{$loc}}{'REF'}\t${$data{$loc}}{'ALT'}\t${$data{$loc}}{'new'}\t${$data{$loc}}{'TRD'}\t${$data{$loc}}{'missing'}[1]\t${$data{$loc}}{'missing'}[0]\t${$data{$loc}}{'allele counts'}[0]\t${$data{$loc}}{'genotype counts'}[0]\n");
    }
    close(OUA);

    print("combined stats file written\n");

} # end of chr loop
