#!/usr/bin/perl

use strict;


my @chrs = (2..19, 'X', 'Y', 'M');
my @mice = (101, 211, 218, 227, 231, 234, 236, 246, 303, 310, 311, 313, 314, 327, 329, 330, 334, 34, 44, 54);


my ($i, $m, $line, $mouse, $ch, $bppos, $ref, $alt, $qch);
my (@row, @data_arr, @geno_arr);
my %snps;


my @indm;
for ($m=0; $m<20; $m++) {
    $indm[$mice[$m]] = $m;
}


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

    my $snpfile = join('', 'SNPs/comb/filter/snps.combined.filter.phr15.chr', $ch, '.vcfstyle.txt');

    my @inds;
    my $Nsnps = 0;
    foreach $mouse ( @mice ) {
#	print("mouse $mouse\n");

	my $fname = join('', 'SNPs/indiv/filter/snps.indiv.filter.phr15.', $mouse, '.chr', $ch, '.sdi.txt');

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

	$line = <IN>;
	chomp($line);
	my $header_line = $line;

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

	    chomp($line);
	    @row = split("\t", $line);
	    $bppos = $row[1];
	    $ref = $row[2];

	    if ( exists( $snps{$bppos} ) ) {

                # check REF alleles match                                                     
                if ( $ref ne ${$snps{$bppos}}{'REF'} ) {
#                   print("ERROR - REF alleles do not match: $ bppos  $ref ${$snps{$bppos}}{'REF'}\n");                     
#                   exit;                                    
                    $snps{$bppos}{'DELETE'} = 'Y';
                }

            } else {
                # add snp to hash and index vector                   
                $snps{ $bppos }{'REF'} = $ref;
		$Nsnps = ( $Nsnps + 1 );
		$inds[$bppos] = $Nsnps;
            }

	    my @tmp = @row[6,8..11];

            if ( $row[4] eq $row[6] ) {
                $qch = 0;
            } else {
                $qch = 1;
            }
	    push(@tmp, $qch);

	    my $tmp_line = join(':', @tmp);
	    $data_arr[$inds[$bppos]][$indm[$mouse]] = $tmp_line;

	} # end of while loop
	close(IN);

    } # end of mouse loop


    my ($pos, $sumc, $alt_str, $gen);
    my (@alts, @alls);


    open(OUS, '>', $snpfile) or die "Cannot open file $snpfile for writing, reason $!\n";

    my $format_header_str1 = '##FORMAT GT=genotype using high mapping quality reads, RD=read depth at variant position, DET=imrdenom detection score, PHR=phred score for snps, QCH=indicator if genotype using all reads differs from GT: 0=no difference, 1=different';
    my $format_header_str2 = '##FORMAT imrdenom detection score: 1 detected by both IMR and DENOM, 2 detected by IMR as heterozygous SNP but detected by DENOM as homozygous SNP, 4 detected by IMR only, 5 detected by DENOM only';
    my @header = ('#CHR', 'POS', 'REF', 'ALT', 'TRD', 'FORMAT', @mice);
    my $header_line = join("\t", @header);

    print(OUS "$format_header_str1\n");
    print(OUS "$format_header_str2\n");
    print(OUS "$header_line\n");

    my $format_str = 'GT:RD:DET:PHR:QCH';

    my $Ntot = 0;
    my $Nalt = 0;
    my $Ndel = 0;
    my $Nnone = 0;

    foreach $pos ( sort { $a <=> $b } keys %snps ) {

	if ( exists( ${$snps{$pos}}{'DELETE'} ) ) {

	    # drop SNP
	    $Ndel = ( $Ndel + 1 );

	} else {

	    $Ntot = ( $Ntot + 1 );

	    my %store;

	    for ($m=0; $m<20; $m++) {

		if ( exists( $data_arr[$inds[$pos]][$m] ) ) {

		    my @mouse_data = split(':', $data_arr[$inds[$pos]][$m]);
		    # 6,8..11,qch ie geno, alt, det, rd, phr, qch

		    # ALTs
		    @alts = split(',', $mouse_data[1]);

		    for ( $i=0; $i<=$#alts; $i++ ) {
			if ( exists( $store{$alts[$i]} ) ) {
			    $store{ $alts[$i] } = ( $store{ $alts[$i] } + 1 );
			} else {
			    $store{ $alts[$i] } = 1;
			}
		    }
		} # end of if loop
	    } # end of sum mice loop


	    my @alt_vec;
	    foreach my $al ( keys %store ) {
		push(@alt_vec, $al);
	    }
	    $alt_str = join(',', @alt_vec);
	    if ( length($alt_str) < 1 ) {
		$alt_str = '.';
	    }
	    $snps{$pos}{'ALT'} = $alt_str;

	    if ( length($alt_str) > 1 ) {
		$Nalt++;
	    }

	    # genotypes

	    my @info_vec;

	    my $Nmiss = 0;
	    my $sumc = 0;

	    for ($m=0; $m<20; $m++) {

		if ( exists( $data_arr[$inds[$pos]][$m] ) ) {

		    my @mouse_data = split(':', $data_arr[$inds[$pos]][$m]);
		    # 6,8..11,qch ie geno, alt, det, rd, phr, qch             

		    my ($al1, $al2, $gen);

		    if ( $mouse_data[0] eq '..' || $mouse_data[0] eq '.' ) {
			$gen = './.';
			$Nmiss++;
		    } else {
			@alls = split('', $mouse_data[0]);

			if ( $alls[0] eq ${$snps{$pos}}{'REF'} ) {
			    $al1 = 0;
			} else {
			    for ($i=0; $i<=$#alt_vec; $i++) {
				if ( $alls[0] eq $alt_vec[$i] ) {
				    $al1 = ($i + 1);
				}
			    }
			} # end of REF loop
			if ( $alls[1] eq ${$snps{$pos}}{'REF'} ) {
			    $al2 = 0;
			} else {

			    for ($i=0; $i<=$#alt_vec; $i++) {
				if ( $alls[1] eq $alt_vec[$i] ) {
				    $al2 = ($i + 1);
				}
			    }
			} # end of REF loop
			if ( defined($al1) && defined($al2) ) {
			    $gen = join('', $al1, '/', $al2);
			} else {
			    print("ERROR - undefined allele: bp $pos mouse $m  mouse data @mouse_data  alleles 1 $al1,   2 $al2  allele vector for genotype $mouse_data[0] is @alls\n");
			    exit;
			}
		    } # end of missing loop

		    # make mouse info string

		    my @mouse_info;
		    push(@mouse_info, $gen);
		    push(@mouse_info, $mouse_data[3]);
		    push(@mouse_info, $mouse_data[2]);
		    push(@mouse_info, $mouse_data[4]);
		    push(@mouse_info, $mouse_data[5]);
		    my $mouse_str = join(':', @mouse_info);
		    push(@info_vec, $mouse_str);


		    # coverage 
                    $sumc = ( $sumc + $mouse_data[3]);

		} else {

		    # no data for this mouse - make mouse info string

		    my @mouse_info;
                    push(@mouse_info, '0/0');
                    push(@mouse_info, '.');
                    push(@mouse_info, '.');
                    push(@mouse_info, '.');
                    push(@mouse_info, '.');
                    my $mouse_str = join(':', @mouse_info);
		    push(@info_vec, $mouse_str);

		} # end of if exists loop

	    } # end of genotypes mouse loop


	    if ( $Nmiss == 20 ) {

		# drop snp - no genotypes

		$Nnone++;

	    } else {

		$snps{$pos}{'TRD'} = 20*$sumc/(20 - $Nmiss);

		# make output line

		my @out_vec = ($ch, $pos, ${$snps{$pos}}{'REF'}, ${$snps{$pos}}{'ALT'}, ${$snps{$pos}}{'TRD'}, $format_str, @info_vec);
		my $out_line = join("\t", @out_vec);
		print(OUS "$out_line\n");
	    }

	} # end of if delete loop
    } # end of foreach pos loop
    close(OUS);

    print("Total snps $Ntot\nDeleted snps (diff REFs) $Ndel\n mult ALTs $Nalt\nsnps with no genotypes $Nnone\n");

    undef %snps;
    undef @data_arr;
} # end of chr loop
