#!/usr/bin/perl

use strict;


my @chrs = (1..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 @mice = (246);
my $phred_thres = 15;


my ($line, $mouse, $ch, $det, $rdep, $phr);
my @row;


# my @header = ('CHR', 'POS', 'REF', 'CALL', 'GT', 'GT.AC', 'GT.MQ', 'GT.MQ.AC', 'ALT', 'DET.SC', 'RDEP', 'PHRED');


foreach $mouse ( @mice ) {
    print("mouse $mouse\n");

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

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

	my $Ntot = 0;
	my $Nrdep = 0;
	my $Ndet = 0;
	my $Nphr = 0;
	my $Ndr = 0;

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

	open(OUS, '>', $ffile) or die "Cannot open file $ffile for writing, reason $!\n";
	my $header_line = <IN>;
	print(OUS "$header_line");

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

	    chomp($line);
	    @row = split("\t", $line);
	    $Ntot = ( $Ntot + 1 );
	    $rdep = $row[10];

	    if ( $rdep < 3 ) {

		$Nrdep++;

	    } else {

		$det = $row[9];
		$phr = $row[11];

		if ( $det == 1 || $det == 2 ) {
		    # keep snp

		    $Ndet++;
		    print(OUS "$line\n");

		} elsif ( $det == 4 || $det == 5 ) {
		    # filter on phred score

		    if ( $phr >= $phred_thres ) {
			# keep snp

			$Nphr++;
			print(OUS "$line\n");

		    } else {
			# drop snp
			$Ndr++;
		    }

		} else {
		    print("ERROR - unrecognised detection score $det\n");
		    exit;
		}
	    }
	} # end of while loop

	close(OUS);
	close(IN);

	print("Total snps $Ntot    Keep det $Ndet    Keep phr $Nphr    Dropped $Ndr    low coverage $Nrdep\n");
    } # end of chr loop
    print("\n\n");
} # end of mouse loop
