
#!/usr/bin/perl

use strict;


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

my @mice = (101, 211, 218, 227, 231, 234, 236, 246, 303, 310, 311, 313, 314, 327, 329, 330, 334, 34, 44, 54);

foreach my $ch ( @chrs ) {

    print("chromosome $ch\n");

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


    my @index = (0) x 200000000;

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

#    open(OUT, '>', $ofile) or die "Cannot open file $ofile for writing, reason $!\n "; 

    my ($line, $start, $bppos, $i, $j);
    my @row;
    my @ind_arr;


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

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

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

    my $Ntot = 0;

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

	$Ntot++;

	for ($i=0; $i<20; $i++) {
	    $ind_arr[$Ntot][$i] = 0;
	}

	@row = split("\t", $line); 
	$bppos = $row[1];
	$index[$bppos] = $Ntot;
    }
    close(INW);


    # read in INDELs files and make indicator array

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

	my $ifile = join('', 'INDELs/wildmus', $mice[$i], '.imrdenom.indels.chr', $ch, '.sdi');

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

	my $Ndel = 0;
	my $Ndsnp = 0;

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

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

	    if ( $row[2] < 0 ) {

		$Ndel++;

		my $first = $row[1];
		my $last = ($row[1] - $row[2] - 1);

		if ( $first == $last ) {

		    if ( $index[$first] > 0 ) {
			# snp is an indel        
			$ind_arr[$index[$j]][$i] = 1;
			$Ndsnp++;
		    }

		} else {

		    for ($j=$first; $j<=$last; $j++) {
			if ( $index[$j] > 0 ) {
			    # snp inside an indel
			    $ind_arr[$index[$j]][$i] = 1;
			    $Ndsnp++;
			}
		    }
		}


	    } # end of deletion check loop
	} # end of while loop
	close(INS);
#	print("mouse $i:  Total snps $Ntot  snps in deletions $Ndsnp  total deletions $Ndel\n");
    } # end of indel read in loop over mice


    # make sum vector

    my @sum_ind = (0) x $Ntot;

    for ($i=0; $i<$Ntot; $i++) {
	for ($j=0; $j<20; $j++) {
	    $sum_ind[$i] = ( $sum_ind[$i] + $ind_arr[$i][$j]);
	}
    }

    open(INW, '<', $fname) or die "Cannot open vcfstyle file $fname, reason $!\n";
    open(OUT, '>', $ofile) or die "Cannot open file $ofile for writing, reason $!\n ";              

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

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

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

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

    my $Nmod = 0;
    my $Nalr = 0;
    my $Ndel = 0;

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

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

	if ( $sum_ind[$index[$bppos]] == 20 ) {

	    # drop snp
	    $Ndel++;

	} elsif ( $sum_ind[$index[$bppos]] > 0 ) {

	    # some mice have an indel at this position
	    my @mouse_genos = @row[6..25];

	    for ($i=0; $i<20; $i++) {
		if ( $ind_arr[$index[$bppos]][$i] == 1 ) {
		    $Nmod++;
		    my @tmp_arr = split(':', $mouse_genos[$i]);

		    if ( $tmp_arr[0] ne './.' ) {
			$tmp_arr[0] = './.';
			my $tmp_str = join(':', @tmp_arr);
			$mouse_genos[$i] = $tmp_str;
		    } else {
			$Nalr++;
		    }
		}
	    }

	    my @new_row = ( @row[0..5], @mouse_genos );
	    my $new_line = join("\t", @new_row);

	    print(OUT "$new_line\n");

	} else {
	    # no snps in indels
	    print(OUT "$line\n");
	}
    }
    close(INW);
    close(OUT);

    undef @ind_arr;
    undef @index;
    undef @sum_ind;

    print("N deletion snps $Ndel  N modified genotypes $Nmod  unnecessary $Nalr\n");

} # end of chr loop
