#!/usr/bin/perl

use strict;



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

my ($line, $ch);

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

    my $INfile = 'MUGAsnps.Flines.TAUlines.chr'.$ch.'.csv';

    my $mapfile = 'GENOTYPES/chr'.$ch.'.muga.map';
    my $allelesfile = 'GENOTYPES/chr'.$ch.'.muga.alleles';
    my $datafile = 'GENOTYPES/chr'.$ch.'.muga.data';


    open(IN, '<', $INfile) or die "Cannot open file $INfile, reason $!\n";
    open(MAP, '>', $mapfile) or die "Cannot open map file $mapfile for writing, reason $!\n";

    my $header_line = <IN>;

    my @header = ('marker', 'chromosome', 'bp');
    my $map_header_line = join("\t", @header);
    print(MAP "$map_header_line\n");


    my ($chr, $map_line, $Nsnps);
    my (@row, @map_row);

    $Nsnps = 0;

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

	$Nsnps = ( $Nsnps + 1 );

	@row = split(',', $line);

	$chr = join('', 'chr', $row[1]);
	@map_row = ($row[0], $chr, $row[2]);
	$map_line = join("\t", @map_row);

	print(MAP "$map_line\n");
    } # end of while loop

    close(MAP);
    close(IN);



    open(ALL, '>', $allelesfile) or die "Cannot open alleles file $allelesfile for writing, reason $!\n";
    print(ALL "markers $Nsnps strains 8\n");

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

    my $header_line = <IN>;
    chomp($header_line);
    my @header = split(',', $header_line);

#    print("@header\n");

    my ($i, $mb);
    my (@strains, @name_vec, @trim);

    for ($i=3; $i<11; $i++) {
	@trim = split(/\*\*\*/, $header[$i]);
	@name_vec = split(/m[0-9]/, $trim[1]);
	push(@strains, $name_vec[0]);
    }
    print(ALL "strain_names @strains\n");

    my @tau_header = @header[11..$#header];
    my $NAstr = "allele NA  0.1250 0.1250 0.1250 0.1250 0.1250 0.1250 0.1250 0.1250";

    my $Nsnps = 0;
    my @geno_arr;

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

	$Nsnps++;

        @row = split(',', $line);
        $chr = join('', 'chr', $row[1]);
	$mb = $row[2]/2000000;
	my @alleles = @row[3..10];
	my @genotypes = @row[11..$#row];

	my %count;
	my $co = 0;
	for ($i=0; $i<8; $i++) {
	    if ( exists($count{$alleles[$i]}) ) {
		$count{$alleles[$i]} = ( $count{$alleles[$i]} + 1 );
	    } else {
		$count{ $alleles[$i]} = 1;
		$co = ( $co + 1 );
	    }
	}
	$co = ( $co + 1 );
	print(ALL "marker $row[0] $co $chr $mb\n");


	foreach my $all ( keys %count ) {

	    if ( $all eq 'A' || $all eq 'C' || $all eq 'G' || $all eq 'T' ) {

		$count{$all} = ( $count{$all}/8 );

		# make output line
		my @prob_vec;
		for ($i=0; $i<=$#alleles; $i++) {
		    if ( $alleles[$i] eq $all ) {
			push(@prob_vec, $count{$all});
		    } else {
			push(@prob_vec, 0);
		    }
		}

		my $prob_str = join(' ', @prob_vec);
		my $out_str = join('  ', "allele $all", $prob_str);
		printf(ALL "allele %s  %.4f %.4f %.4f %.4f %.4f %.4f %.4f %.4f\n", $all, $prob_vec[0], $prob_vec[1], $prob_vec[2], $prob_vec[3], $prob_vec[4], $prob_vec[5], $prob_vec[6], $prob_vec[7]);
	    } else {
		print("ERROR unrecognised founder strain allele $all\n");
		exit;
	    }
	}
	print(ALL "$NAstr\n");


	my @snp_alls = keys %count;
	if ( $#snp_alls > 1 ) {
	    exit;
	}

	for ($i=0; $i<=$#genotypes; $i++) {
	    if ( exists($count{$genotypes[$i]}) ) {
		$geno_arr[$i][(2*$Nsnps - 2)] = $genotypes[$i];
		$geno_arr[$i][(2*$Nsnps - 1)] = $genotypes[$i];
	    } elsif ( $genotypes[$i] eq 'H' ) {
		$geno_arr[$i][(2*$Nsnps - 2)] = $snp_alls[0];
                $geno_arr[$i][(2*$Nsnps - 1)] = $snp_alls[1];
	    } elsif ( $genotypes[$i] eq 'N' ) {
		$geno_arr[$i][(2*$Nsnps - 2)] = 'NA';
                $geno_arr[$i][(2*$Nsnps - 1)] = 'NA';
	    } else {
		print("genotype $genotypes[$i]\n");
		exit;
	    }
	}
    }
    close(IN);
    close(ALL);

    print("N snps $Nsnps\n");


    if ( $Nsnps > 0 ) {

	open(DAT, '>', $datafile) or die "Cannot open data file $datafile for writing, reason $!\n";

	for ($i=0; $i<=$#tau_header; $i++) {

	    my @tmp = split(/\*\*\*/, $tau_header[$i]);

	    if ( $tmp[1] =~ m/IL/ ) {
		if ( $tmp[1] =~ m/m/ ) {
		    my @ilstr = split('m', $tmp[1]);
		    my $ilstring = join('', substr($ilstr[0], 0, 2), '-', substr($ilstr[0], 2));

		    my $out_str = join(' ', "IL", $ilstring, "0 0 M NA", @{$geno_arr[$i]});
		    print(DAT "$out_str\n");

		} elsif ( $tmp[1] =~ m/f/ ) {
		    my @ilstr = split('f', $tmp[1]);
		    my $ilstring = join('', substr($ilstr[0], 0, 2), '-', substr($ilstr[0], 2));

		    my $out_str = join(' ', "IL", $ilstring, "0 0 F NA", @{$geno_arr[$i]});
		    print(DAT "$out_str\n");

		} else {
		    print("no m or f in 2nd half\n");
		    exit;
		}
	    } elsif ( $tmp[0] =~ m/IL/ ) {

		if ( $tmp[1] =~ m/m/ ) {
		    my $ilstring = join('', substr($tmp[0], 0, 2), '-', substr($tmp[0], 2));

		    my $out_str = join(' ', "IL", $ilstring, "0 0 M NA", @{$geno_arr[$i]});
		    print(DAT "$out_str\n");

		} elsif ( $tmp[1] =~ m/f/ ) {
		    my $ilstring = join('', substr($tmp[0], 0, 2), '-', substr($tmp[0], 2));

		    my $out_str = join(' ', "IL", $ilstring, "0 0 F NA", @{$geno_arr[$i]});
		    print(DAT "$out_str\n");

		} else {
		    print("no m or f in 2nd half\n");
		    exit;
		}
	    } else {
		print("no IL string @tmp\n");
		exit;
	    }
	}
	close(DAT);

    } # end of if Nsnps loop

} # end of chr loop

