
#!/usr/bin/perl

use strict;


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


my $snpset = $ARGV[0];
my $subdir = $ARGV[1];


my $sumfile = $subdir.'/crossspecies.sumout.txt';

open(OUT, '>', $sumfile) or die "Cannot open file $sumfile for writing, reason $!\n";
print(OUT "chr N.dom N.cast N.both N.dom.uniq N.cast.uniq N.m.ref N.m.ref.check\n");


foreach my $ch ( @chrs ) {

    print("chromosome $ch\n");

    my $dom_name = '/Net/sparse/data/wildmice/callsets/domGATK'.$snpset.'/snps.filtered.genotypes.filtered.chr'.$ch.'.sub.dom.vcf';
    my $cast_name = '/Net/sparse/data/wildmice/callsets/castGATK'.$snpset.'/snps.filtered.genotypes.filtered.chr'.$ch.'.sub.cast.vcf';

    my $both_dom_file = $subdir.'/SNPs.both.dom.chr'.$ch.'.vcf';
    my $both_cast_file = $subdir.'/SNPs.both.cast.chr'.$ch.'.vcf';
    my $uniq_dom_file = $subdir.'/SNPs.uniq.dom.chr'.$ch.'.vcf';
    my $uniq_cast_file = $subdir.'/SNPs.uniq.cast.chr'.$ch.'.vcf';


    my ($line, $start, $bppos, $ref);
    my @row;
    my (%dom, %cast);


    open(IND, '<', $dom_name) or die "Cannot open vcf file $dom_name, reason $!\n";

    $line = <IND>;
    $start = substr( $line, 0, 2);

    while ( $start eq '##' ) {
        $line = <IND>;
        $start = substr( $line, 0, 2);
    } # end of headers while loop                                                                                     
 
    my $header_line = $line;
    chomp($header_line);

    my $Ndom = 0;

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

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

	if ( exists( $dom{$bppos} ) ) {
	    print("location already in dom hash $bppos\n");
	} else {
	    $dom{ $bppos }{'REF'} = $ref;
	    $Ndom++;
	}
    }
    close(IND);
    print("hash of dom snps created\n");


    open(OUBC, '>', $both_cast_file) or die "Cannot open file $both_cast_file for writing, reason $!\n";
    open(OUOC, '>', $uniq_cast_file) or die "Cannot open file $uniq_cast_file for writing, reason $!\n";
    open(INC, '<', $cast_name) or die "Cannot open vcf file $cast_name, reason $!\n";

    $line = <INC>;
    my $start = substr( $line, 0, 2);
    print(OUBC "$line");
    print(OUOC "$line");

    while ( $start eq '##' ) {
	$line = <INC>;
	$start = substr( $line, 0, 2);
	print(OUBC "$line");
	print(OUOC "$line");
    } # end of headers while loop  

    my $header_line = $line;
    chomp($header_line);

    my $Ncast = 0;
    my $Nonlycast = 0;
    my $Nboth = 0;
    my $Nmref = 0;

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

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

	if ( exists( $cast{$bppos} ) ) {
	    print("location already in cast hash $bppos\n");
	} else {

	    $cast{ $bppos }{'REF'} = $ref;
	    $Ncast++;

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

		# check refs match

		if ( ${$dom{$bppos}}{'REF'} eq ${$cast{$bppos}}{'REF'} ) {

		    $Nboth++;
		    print(OUBC "$line\n");

		} else {
		    $Nmref++;
		} # end of ref match loop
	    } else {
		# SNP unique to cast
		$Nonlycast++;
		print(OUOC "$line\n");
	    } # end of if both loop

	} # end of if duplicate lines loop
    } # end of while loop
    close(INC);
    close(OUBC);
    close(OUOC);

    print("Cast file read and hash of cast snps created\n");


    open(OUBD, '>', $both_dom_file) or die "Cannot open file $both_dom_file for writing, reason $!\n ";
    open(OUOD, '>', $uniq_dom_file) or die "Cannot open file $uniq_dom_file for writing, reason $!\n ";


    open(IND, '<', $dom_name) or die "Cannot open vcf file $dom_name, reason $!\n";

    $line = <IND>;
    $start = substr( $line, 0, 2);

    while ( $start eq '##' ) {
	print(OUBD "$line");
	print(OUOD "$line");

	$line = <IND>;
        $start = substr( $line, 0, 2);
    } # end of headers while loop                                                                                     

    my $header_line = $line;
    print(OUBD "$header_line");
    print(OUOD "$header_line");

    my $Nboth = 0;
    my $Nonlydom = 0;
    my $Nmref_cast = 0;

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

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

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

	    # check refs match

	    if ( $ref eq ${$cast{$bppos}}{'REF'} ) {

		$Nboth++;
		print(OUBD "$line\n");

	    } else {
		$Nmref_cast++;
	    }

	} else {
	    $Nonlydom++;
	    print(OUOD "$line\n");
	}
    }
    close(IND);

    close(OUBD);
    close(OUOD);

    print(OUT "$ch $Ndom $Ncast $Nboth $Nonlydom $Nonlycast $Nmref $Nmref_cast\n");


#    print("number of snps $Ndom\n");
#    print("number of snps $Ncast\nOverlap snps $Nboth, snps unique to cast $Nonlycast\nsnps with REF mismatches $Nmref\n");
#    print("Dom file read: overlap snps $Nboth should match cast above, snps in dom but not cast $Nonlydom\nsnps with REF mismatches with cast $Nmref_cast\n");

} # end of chr loop

close(OUT);
