#!/usr/bin/perl

use strict;


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

my @prel_indices = (0, 1, 3, 4);
#my @data_indices = (9..28);
my @data_indices = (9,12,13,18,19,24,25,10,11,15,14,16,17,21,20,22,23,27,26,28);


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('', 'vcf/SNPs.known.chr', $ch, '.vcf');
    print("reading $fname\n");

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

    open(OUT, '>', 'alleles/alleles.known.chr'.$ch.'.txt') or die "Cannot open file 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[3];
	$alt = $row[4];

	@mice = @row[ @data_indices ];

	my @alleles;
	my @depths;

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

	    $qual = $mouse_info[2];

	    @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 $fname = join('', 'vcf/SNPs.new.chr', $ch, '.vcf');
    print("reading $fname\n");

    my %new;

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

    open(OUT, '>', 'alleles/alleles.new.chr'.$ch.'.txt') or die "Cannot open file 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[3];
	$alt = $row[4];

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

	@mice = @row[ @data_indices ];

	my @alleles;
	my @depths;

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

	    $qual = $mouse_info[2];

	    @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 $fname = join('', 'vcf/snps.filtered.hwe.chr', $ch, '.vcf');
    print("reading $fname\n");

    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], 'NEW', @GT, @RD);
#    print("$new_header_line\n");

    open(OUT, '>', 'alleles/alleles.chr'.$ch.'.txt') or die "Cannot open file 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[3];
	$alt = $row[4];

	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);
	    
	    $qual = $mouse_info[2];

	    @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);

} # end of chr loop
