
#!/usr/bin/perl -w

use strict;


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

foreach my $ch ( @chrs ) {

    print("\nchromosome $ch\n");

    my $gatk_fname = join('', '../gatk/alleles/alleles.chr', $ch,'.txt');
    my $plat_fname = join('', '../platypus/alleles/alleles.chr', $ch, '.txt');

    print("gatk input file $gatk_fname\n");
    print("plat input file $plat_fname\n");

    my $bothfile = join('', '../intersection.data/gatk.plat/alleles.both.chr', $ch, '.txt');
    my $gatkfile = join('', '../intersection.data/gatk.plat/alleles.onlygatk.chr', $ch, '.txt');
    my $platfile = join('', '../intersection.data/gatk.plat/alleles.onlyplat.chr', $ch, '.txt');


    my ($line, $bppos);
    my (@row);
    my (%gatk, %plat, %both);


    open(OUTB, '>', $bothfile) or die "Cannot open file $bothfile for writing, reason $!\n ";
    open(INW, '<', $gatk_fname) or die "Cannot open alleles file $gatk_fname, reason $!\n";

    $line = <INW>;
    my $gatk_header_line = $line;

    my $Ngatk = 0;

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

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

	if ( exists( ${$gatk{$bppos}}{'data'} ) ) {
	    print("location already in gatk hash  $bppos\n");
	    exit;
	} else {
	    $gatk{ $bppos }{'data'} = $line;
	    $Ngatk = ( $Ngatk + 1 );
	}
    }
    close(INW);
    print("hash of gatk data created\nNumber of SNPs\n");



    my $i;
    my ($gatk_ref, $plat_ref, $gatk_alt, $plat_alt, $ALT_match);

    open(OUTP, '>', $platfile) or die "Cannot open file $platfile for writing, reason $!\n ";
    open(INS, '<', $plat_fname) or die "Cannot open alleles file $plat_fname, reason $!\n";

    $line = <INS>;
    my $plat_header_line = $line;
    print(OUTP "$plat_header_line");

    # make output header for both file

    chomp($plat_header_line);
    chomp($gatk_header_line);
    my @plat_header = split("\t", $plat_header_line);
    my @gatk_header = split("\t", $gatk_header_line);
    my @out_header = @plat_header[0..2];

    $out_header[3] = "gatk.$gatk_header[3]";
    $out_header[4] = "plat.$plat_header[3]";
    $out_header[5] = 'ALT.match';
    $out_header[6] = "gatk.$gatk_header[4]";
    $out_header[7] = "plat.$plat_header[4]";
    $out_header[8] = 'Num.geno.diffs';

    for ($i=5; $i<25; $i++) {
	$out_header[(4 + $i)] = join('.', 'match', $gatk_header[$i]);
	$out_header[(24 + $i)] = join('.', 'gatk', $gatk_header[$i]);
	$out_header[(44 + $i)] = join('.', 'plat', $plat_header[$i]);
    }
    my $out_header_line = join("\t", @out_header);
    print(OUTB "$out_header_line\n");


    my $Nplat = 0;
    my $Nonlyplat = 0;
    my $Nboth = 0;
    my $Naltwr = 0;
    my $Nrefri = 0;
    my $Nmatch = 0;

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

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

	if ( exists( ${$plat{$bppos}}{'data'} ) ) {
	    print("location already in plat hash  $bppos\n");
	    exit;
	} else {

	    $plat{ $bppos }{'data'} = $line;
	    $Nplat = ( $Nplat + 1 );

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

		$Nboth = ( $Nboth + 1 );

		# check SNP data matches
		my @gatk_vec = split("\t", ${$gatk{$bppos}}{'data'} );
		$gatk_ref = $gatk_vec[2];
		$gatk_alt = $gatk_vec[3];

		my @plat_vec = split("\t", ${$plat{$bppos}}{'data'} );
		$plat_ref = $plat_vec[2];
		$plat_alt = $plat_vec[3];

		if ( $gatk_ref ne $plat_ref ) {
		    print("ERROR - REFs do not match  pos = $bppos  REFs gatk $gatk_ref  plat $plat_ref\n");
		    exit;
		} else {
		    # SNP called by both callers 

		    $Nrefri = ( $Nrefri + 1 );

		    $gatk{$bppos}{'both'} = 'Y';

		    # check ALTs match

		    if ( $gatk_alt ne $plat_alt ) {
#			print("ALTs do not match  pos = $bppos  ALTs gatk $gatk_alt  plat $plat_alt\n");
			$Naltwr = ( $Naltwr + 1 );
			$ALT_match = 'N';
		    } else {
			$ALT_match = 'Y';
		    }

		    # check genotypes match for each mouse

		    my @gatk_genos = @gatk_vec[5..24];
		    my @plat_genos = @plat_vec[5..24];
		    my @matches;

		    my $Nmism = 0;

		    for ($i=0; $i<20; $i++) {
#			print("mouse $i: $gatk_genos[$i] $plat_genos[$i]\n");
			if ( $gatk_genos[$i] eq '.' || $plat_genos[$i] eq '.') {
			    push( @matches, '.');
			} elsif ( $gatk_genos[$i] == $plat_genos[$i] ) {
			    push( @matches, 'Y');
			} else {
			    push( @matches, 'N');
			    $Nmism = ( $Nmism + 1 );
			}
		    }
#		    print("Number of mismatches = $Nmism\n");
		    if ( $Nmism == 0 ) {
			$Nmatch = ( $Nmatch + 1 );
		    }

		    # make output vector
		    # gatk_vec, plat_vec

		    my @out_vec = @gatk_vec[0..3];
		    $out_vec[4] = $plat_vec[3];
		    $out_vec[5] = $ALT_match;
		    $out_vec[6] = $gatk_vec[4];
		    $out_vec[7] = $plat_vec[4];
		    $out_vec[8] = $Nmism;
		    for ($i=0; $i<20; $i++) {
			$out_vec[(9 + $i)] = $matches[$i];
			$out_vec[(29 + $i)] = $gatk_vec[(5 + $i)];
			$out_vec[(49 + $i)] = $plat_vec[(5 + $i)];
		    }
		    my $out_line = join("\t", @out_vec);

		    print(OUTB "$out_line\n");

		} # ref match loop

	    } else {
		# SNP only called by platypus
		print(OUTP "$line\n");
		$Nonlyplat = ( $Nonlyplat + 1 );
	    } # gatk loop
	} # check duplicates loop
    } # while loop
    close(INS);
    close(OUTP);
    close(OUTB);
    print("hash of plat data created\noutput only PLAT file written\n");

    print("Number of SNPs GATK $Ngatk  PLAT $Nplat\n Number only in PLAT $Nonlyplat\n Number in both $Nboth\n Number with ref matches $Nrefri\n Number with ALT mismatches $Naltwr\nNumber of SNPs in both with all genotype matches $Nmatch\n");


    open(OUTG, '>', $gatkfile) or die "Cannot open file $gatkfile for writing, reason $!\n ";
    print(OUTG "$gatk_header_line\n");

    my $loc;

    my @gatk_srted_keys;
    @gatk_srted_keys = sort { $a <=> $b } keys( %gatk );

    my $Cboth = 0;
    my $Nonlygatk = 0;
    my $Nerrors = 0;

    foreach $loc ( @gatk_srted_keys ) {

	if ( exists ( ${$gatk{ $loc }}{'both'} ) ) {
	    # skip SNP - in both data sets
	    $Cboth = ( $Cboth + 1 );
	} else {
	    if ( exists( ${$gatk{ $loc }}{'data'} ) ) {
		$Nonlygatk = ( $Nonlygatk + 1 );
		print(OUTG "${$gatk{ $loc }}{'data'}\n");
	    } else {
		$Nerrors = ( $Nerrors + 1 );
	    }
	}
    }

    close(OUTG);
    print("output only GATK file written\n");

    print("Number of SNPs only in GATK $Nonlygatk\nNumber of error keys $Nerrors\nCheck number in both $Cboth\n");

} # end of chr loop
