
#!/usr/bin/perl

use strict;


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

foreach my $ch ( @chrs ) {

    print("chromosome $ch\n");

    my $wm_fname = join('', 'vcf/snps.filtered.hwe.chr', $ch,'.vcf');
    my $s17_fname = join('', '../../../mouse.snps/17strain.data/final.chr', $ch, '.vcf');

    print("input file $wm_fname\n");
    print("17 strains file $s17_fname\n");

    my $knownfile = join('', 'vcf/SNPs.known.chr', $ch, '.vcf');
    my $newfile = join('', 'vcf/SNPs.new.chr', $ch, '.vcf');
    my $known17strfile = join('', 'vcf/SNPs.known.17straindata.chr', $ch, '.vcf');


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

    my ($line, $start);

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

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

    my ($line, $bppos, $ref, $alt, $mice_string);
    my (@row, @mice);
    my %wild;

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

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

	$bppos = $row[1];

	if ( exists( ${$wild{$bppos}}{'data'} ) ) {
	    print("location already in hash  $bppos\n");
	} else {
	    $wild{ $bppos }{'data'} = 'Y';
	}
    }
    close(INW);

    print("hash of wild mouse data created\n");


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

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

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

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

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

    my $header_line = $line;
    print(OUTS "$header_line");


    my $N17 = 0;

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

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

	if ( exists( ${$wild{ $bppos }}{'data'} ) ) {

	    $N17 = ( $N17 + 1 );

	    print(OUTS "$line");
	    $wild{ $bppos }{'known'} = 'Y';
	}

    } # end of while loop
    close(INS);
    close(OUTS);

    print("matches in 17 strain data found and printed  $N17\n");


    open(OUTK, '>', $knownfile) or die "Cannot open file $knownfile for writing, reason $!\n ";
    open(OUTN, '>', $newfile) or die "Cannot open file $newfile for writing, reason $!\n ";

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

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

    while ( $start eq '##' ) {
	print(OUTK "$line");
	print(OUTN "$line");

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

    my $header_line = $line;
    print(OUTK "$header_line");
    print(OUTN "$header_line");


    my $Nknown = 0;
    my $Nnew = 0;

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

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

	if ( exists( ${$wild{ $bppos }}{ 'known' }) ) {
	    $Nknown = ( $Nknown + 1 );
	    print(OUTK "$line");
	} else {
	    $Nnew = ( $Nnew + 1 );
	    print(OUTN "$line");
	}
    }
    close(INW);

    close(OUTK);
    close(OUTN);

    print("Chromosome $ch: counts of SNPs:\n17 strains $N17\nknown $Nknown\nnew $Nnew\n");

} # end of chr loop
