#!/usr/bin/perl

use strict;


my @chrs = (1..19, 'X', 'Y');
#my @chrs = (1..19);
#print("@chrs\n");
my @indices = (0, 1, 3, 4, 9..25);


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

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


print("reading 17 strains data\n\n");

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

    my $fname = join('', 'vcf/SNPs.known.17straindata.chr', $ch, '.vcf');


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

    my @header = split("\t", $header_line);
    my @new_header = @header[@indices];
    my $new_header_line = join("\t", @new_header);

    open(OUT, '>', 'alleles/alleles.17strains.chr'.$ch.'.txt') or die "Cannot open file for writing, reason $!\n";
    print(OUT "$new_header_line\n");


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

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

	$alt = $row[4];

	if ( length($alt) > 1 ) {

#	    print("skip this SNP: too many alt alleles $alt\n");
	    $nskip = ( $nskip + 1 );

	} else {

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

	    @mice = @row[ 9..$#row ];

	    my @alleles;

	    foreach $mouse ( @mice ) {
		@mouse_info = split(/:/, $mouse);
		@geno_vec = split(/\//, $mouse_info[0]);
		$qual = $mouse_info[1];

		if ( $qual >= 0 ) {
		    $geno = $qual;
		} elsif ( $qual < 0 && $qual > -7 ) {
		    if ( $geno_vec[0] == $geno_vec[1] ) {
			$geno = $geno_vec[0];
		    } else {
			print("ERROR - unrecognised genotype @geno_vec and qual score $qual pair\n");
			exit;
		    }
		} elsif ( $qual == -7 || $qual == -8 ) {
		    $geno = '.';
		} else {
		    print("ERROR - unrecognised genotype @geno_vec with qual score $qual at pos $bppos\n");
		    exit;
		}

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

	    $allele_string = join("\t", @alleles);
	    print(OUT "$chr\t$bppos\t$ref\t$alt\t$allele_string\n");
	} # end of skip SNP loop
    } # end of while loop
    close(IN);
    close(OUT);

    print("skipped $nskip SNPs due to multiple ALT alleles\n");

} # end of chr loop
