#!/usr/bin/perl

use strict;


my @chrs = (1..19, 'X', 'Y');
#print("@chrs\n");
my @prel_indices = (0..3);
my @data_indices = (6..25);


my ($fname, $line, $start, $ch, $nskip);
my ($chr, $bppos, $ref, $alt, $mouse, $qual, $geno, $allele_string, $depth_string);
my ($gt_id, $rd_id, $new_ind);

my (@row, @mice, @mouse_info, @geno_vec);



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

    my $fname = join('', 'vcfstyle/genotypes.known.chr', $ch, '.vcfstyle.txt');
    my $nname = join('', 'vcfstyle/genotypes.new.chr', $ch, '.vcfstyle.txt');
    my $aname = join('', 'vcfstyle/snps.all.biALTs.chr', $ch, '.vcfstyle.txt');

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


    open(IN, '<', $fname) or die "Cannot open file $fname, 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);
#    print("$header_line\n");
    my @header = split("\t", $header_line);

    my @GT;
    my @RD;
    foreach my $id ( @header[@data_indices] ) {
	$gt_id = join('.', 'GT', $id);
	push(@GT, $gt_id);
	$rd_id = join('.', 'RD', $id);
	push(@RD, $rd_id);
    }
    my $new_header_line = join("\t", @header[@prel_indices], @GT, @RD);
#    print("$new_header_line\n");

    my $Nmult_knw = 0;
    my $Nmiss_ALT_knw = 0;

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

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

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

	$chr = $row[0];
	$bppos = $row[1];
	$ref = $row[2];
	$alt = $row[3];

	if ( length($alt) < 1 ) {
	    print("ALT $alt missing\n");
	    $alt = '.';
	    $Nmiss_ALT_knw++;
	}

	if ( length($alt) > 1 ) {

	    $Nmult_knw++;

	} else {

	    @mice = @row[ @data_indices ];

	    my @alleles;
	    my @depths;

	    foreach $mouse ( @mice ) {
		@mouse_info = split(/:/, $mouse);

		if ( exists($mouse_info[1]) ) {
		    $qual = $mouse_info[1];
		} else {
		    $qual = '.';
		}

		@geno_vec = split(/\//, $mouse_info[0]);
		if ( $geno_vec[0] eq '.' ) {
		    if ( $geno_vec[1] eq '.' ) {
			$geno = '.';
			$qual = '.';
		    } else {
			print("missing allele @geno_vec\n");
			exit;
		    }
		} elsif ( $geno_vec[0] == 0 || $geno_vec[0] == 1 ) {
		    if ( $geno_vec[1] eq '.' ) {
			print("missing allele @geno_vec\n");
			exit;
		    } elsif ( $geno_vec[1] == 0 || $geno_vec[1] == 1 ) {
			$geno = ( $geno_vec[0] + $geno_vec[1] );
		    } else {
			print("unrecognised allele $geno_vec[1]\n");
			exit;
		    }
		} else {
		    print("unrecognised allele $geno_vec[0]\n");
		    exit;
		}

		push(@alleles, $geno);
		push(@depths, $qual);
	    } # end of loop over mice

	    $allele_string = join("\t", @alleles);
	    $depth_string = join("\t", @depths);

	    print(OUT "$chr\t$bppos\t$ref\t$alt\t$allele_string\t$depth_string\n");

	}

    } # end of while loop

    close(IN);
    close(OUT);



    my %new;

    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);
    } # end of headers while loop

    my $header_line = substr( $line, 1);
    chomp($header_line);
#    print("$header_line\n");
    my @header = split("\t", $header_line);

    my @GT;
    my @RD;
    foreach my $id ( @header[@data_indices] ) {
	$gt_id = join('.', 'GT', $id);
	push(@GT, $gt_id);
	$rd_id = join('.', 'RD', $id);
	push(@RD, $rd_id);
    }
    my $new_header_line = join("\t", @header[@prel_indices], @GT, @RD);
#    print("$new_header_line\n");

    my $Nmult_new = 0;
    my $Nmiss_ALT_new = 0;

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

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

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

	$chr = $row[0];
	$bppos = $row[1];
	$ref = $row[2];
	$alt = $row[3];

	if ( length($alt) < 1 ) {
            print("ALT $alt missing\n");
	    $alt = '.';
            $Nmiss_ALT_new++;
        }

	if ( length($alt) > 1 ) {

	    $Nmult_new++;

	} else {

	    $new{ $bppos }{'REF'} = $ref;

	    @mice = @row[ @data_indices ];

	    my @alleles;
	    my @depths;

	    foreach $mouse ( @mice ) {
		@mouse_info = split(/:/, $mouse);

		if ( exists($mouse_info[1]) ) {
		    $qual = $mouse_info[1];
		} else {
		    $qual = '.';
		}

		@geno_vec = split(/\//, $mouse_info[0]);
		if ( $geno_vec[0] eq '.' ) {
		    if ( $geno_vec[1] eq '.' ) {
			$geno = '.';
			$qual = '.';
		    } else {
			print("missing allele @geno_vec\n");
			exit;
		    }
		} elsif ( $geno_vec[0] == 0 || $geno_vec[0] == 1 ) {
		    if ( $geno_vec[1] eq '.' ) {
			print("missing allele @geno_vec\n");
			exit;
		    } elsif ( $geno_vec[1] == 0 || $geno_vec[1] == 1 ) {
			$geno = ( $geno_vec[0] + $geno_vec[1] );
		    } else {
			print("unrecognised allele $geno_vec[1]\n");
			exit;
		    }
		} else {
		    print("unrecognised allele $geno_vec[0]\n");
		    exit;
		}

		push(@alleles, $geno);
		push(@depths, $qual);
	    } # end of loop over mice

	    $allele_string = join("\t", @alleles);
	    $depth_string = join("\t", @depths);

	    print(OUT "$chr\t$bppos\t$ref\t$alt\t$allele_string\t$depth_string\n");

	}

    } # end of while loop

    close(IN);
    close(OUT);


    open(IN, '<', $aname) or die "Cannot open file $aname, 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);
#    print("$header_line\n");
    my @header = split("\t", $header_line);

    my @GT;
    my @RD;
    foreach my $id ( @header[@data_indices] ) {
	$gt_id = join('.', 'GT', $id);
	push(@GT, $gt_id);
	$rd_id = join('.', 'RD', $id);
	push(@RD, $rd_id);
    }
    my $new_header_line = join("\t", @header[@prel_indices], 'NEW', @GT, @RD);
#    print("$new_header_line\n");

    my $Nmult_all = 0;
    my $Nmiss_ALT_all = 0;

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

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

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

	$chr = $row[0];
	$bppos = $row[1];
	$ref = $row[2];
	$alt = $row[3];

	if ( length($alt) < 1 ) {
            print("ALT $alt missing\n");
	    $alt = '.';
            $Nmiss_ALT_all++;
        }


	if ( length($alt) > 1 ) {

	    $Nmult_all++;

	} else {

	    if ( exists( $new{ $bppos } )) {
		$new_ind = 'Y';
	    } else {
		$new_ind = 'N';
	    }

	    @mice = @row[ @data_indices ];

	    my @alleles;
	    my @depths;

	    foreach $mouse ( @mice ) {
		@mouse_info = split(/:/, $mouse);

		if ( exists($mouse_info[1]) ) {
		    $qual = $mouse_info[1];
		} else {
		    $qual = '.';
		}

		@geno_vec = split(/\//, $mouse_info[0]);
		if ( $geno_vec[0] eq '.' ) {
		    if ( $geno_vec[1] eq '.' ) {
			$geno = '.';
			$qual = '.';
		    } else {
			print("missing allele @geno_vec\n");
			exit;
		    }
		} elsif ( $geno_vec[0] == 0 || $geno_vec[0] == 1 ) {
		    if ( $geno_vec[1] eq '.' ) {
			print("missing allele @geno_vec\n");
			exit;
		    } elsif ( $geno_vec[1] == 0 || $geno_vec[1] == 1 ) {
			$geno = ( $geno_vec[0] + $geno_vec[1] );
		    } else {
			print("unrecognised allele $geno_vec[1]\n");
			exit;
		    }
		} else {
		    print("unrecognised allele $geno_vec[0]\n");
		    exit;
		}

		push(@alleles, $geno);
		push(@depths, $qual);
	    } # end of loop over mice

	    $allele_string = join("\t", @alleles);
	    $depth_string = join("\t", @depths);

	    print(OUT "$chr\t$bppos\t$ref\t$alt\t$new_ind\t$allele_string\t$depth_string\n");

	}

    } # end of while loop

    close(IN);
    close(OUT);

    print("N mult ALTs: known $Nmult_knw  new $Nmult_new   all $Nmult_all\n");
    print("N missing ALTs: known $Nmiss_ALT_knw  new $Nmiss_ALT_new  all $Nmiss_ALT_all\n");

} # end of chr loop
