#!/usr/bin/perl

use strict;


my $output_selection = 'CONCISE'; # CONCISE or VERBOSE
my $marker_selection = 'PRIV'; # PRIV or CALL


if ($marker_selection eq 'PRIV') {
    print("marker selection = PRIV: uncertain ATG scores treated as PRIVate alleles\n");
} elsif ($marker_selection eq 'CALL') {
    print("marker selection = CALL: uncertain ATG scores are used to CALL probable alleles\n");
} else {
    print("marker selection $marker_selection not recognised: choose either PRIV or CALL\n");
    exit;
}


my @chrs_CC = 1..19;
push (@chrs_CC, 'X');

my @CC_indices = (8, 11, 5, 16, 17, 12, 18, 20);

foreach my $chr ( @chrs_CC ) {

    my $fname = join('', '/data/mouse.snps/17strain.data/final.chr', $chr, '.vcf.txt');
    print("\n$fname\n");

    my $outfile;
    if ($marker_selection eq 'PRIV') {
	$outfile = join('', './PRIV/chr', $chr, '.CC.txt');
	print("output file $outfile\n");
    } elsif ($marker_selection eq 'CALL') {
	$outfile = join('', './CALL/chr', $chr, '.CC.txt');
        print("output file $outfile\n");
    }


    my $i;
    my $Ntotal = 0;
    my $NCC = 0;

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

    my $line = <INF1>;
    $line = <INF1>;
    chomp($line);
    my @header = split(/\t/, $line);

    my @CC_strains = @header[ @CC_indices ];
    my $out_header = join("\t", 'chr', 'bp', @CC_strains);


    open(OUTF, '>', $outfile) or die "Cannot open file $outfile for writing, reason $!\n";
    print( OUTF "$out_header\n");


    while ( defined(my $line = <INF1>) ) {

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

	$Ntotal = ( $Ntotal + 1 );

#  get alleles for CC strains

	my @input_codes;
	foreach my $index ( @CC_indices ) {
	    push(@input_codes, $row[$index]);
	}

	if ( $output_selection eq 'VERBOSE' ) {
            print("alleles in the CC: @input_codes\n");
	}

#  is the marker polymorphic in the CC strains?

	my $same = 0;
	my $Nstr = scalar @input_codes;

	for ($i=0; $i<=$#input_codes; $i++) {
	    if ( $input_codes[$i] == $input_codes[0] ) {
		$same = ($same + 1);
	    }
	}

	if ( $output_selection eq 'VERBOSE' ) {
	    print("same = $same, compare to $Nstr\n");
	}

	if ( $same == $Nstr ) {

            if ( $output_selection eq 'VERBOSE' ) {
                print("Not polymorphic in the CC: skip this marker\n\n");
            }

        } else {

#  get alt allele(s) - are there more than 1?

            my @alt = split(',', $row[3]);
            my $Nalt = scalar @alt;

	    if ( $output_selection eq 'VERBOSE' ) {
		print("N alt alleles @alt = $Nalt\n");
	    }

# process @input_codes

	    my @all_codes;
	    my @qual_codes;
	    my $neg_count = 0;
	    my $unc_calls = 0;
	    foreach my $code ( @input_codes ) {

		my @code_vec = split(/:/, $code);

		push( @all_codes, $code_vec[0] );
		push( @qual_codes, $code_vec[1] );

		if ( $code_vec[1] < 0 ) {
		    $neg_count = ($neg_count + 1);
		    if ( $code_vec[1] < -1 && $code_vec[1] > -8 ) {
			$unc_calls = ($unc_calls + 1);
		    }
		}
	    }

	    if ( $output_selection eq 'VERBOSE' ) {
                print("allele codes @all_codes\n");
		print("quality codes @qual_codes\n");
		print("neg count = $neg_count\n");
		print("unc calls = $unc_calls\n");
	    }

	    if ( $marker_selection eq 'CALL' && $Nalt > 1 && $unc_calls > 0 ) {

		if ( $output_selection eq 'VERBOSE' ) {
		    print("Multi-allelic SNP: ALT alleles cannot be called: skip this marker\n\n");
		}

	    } else {

# substitute negative quality scores

		if ( $neg_count > 0 ) {

		    if ( $marker_selection eq 'PRIV' ) {

			for ($i=0; $i<= $#qual_codes; $i++) {
			    if ( $qual_codes[$i] < 0 ) {
				$all_codes[$i] = -8;
			    }
			}

		    } elsif ( $marker_selection eq 'CALL' ) {

			my $alt_allele;

			for ($i=0; $i<= $#qual_codes; $i++) {
			    if ( $qual_codes[$i] == -1 ) {
				# call as ref
				$all_codes[$i] = 0;
			    } elsif ( $qual_codes[$i] == -8 ) {
				# leave as neg score
				$all_codes[$i] = -8;
			    } elsif ( $qual_codes[$i] == -7 ) {
				# call as het - leave as neg
			    } elsif ( $all_codes[$i] < -1 && $all_codes[$i] > -7 ) {
				# scores between -2 and -6 called as alt
			    }
			}
		    }
		} # end of neg count if loop

		if ( $output_selection eq 'VERBOSE' ) {
		    print("substituted alleles: @all_codes\n");
		}

#  is the marker still polymorphic in the CC strains?

		$neg_count = 0;
		$same = 0;
		$Nstr = scalar @all_codes;

		for ($i=0; $i<=$#all_codes; $i++) {
		    if ( $all_codes[$i] < 0 ) {
			$neg_count = ($neg_count + 1);
		    }
		    if ( $all_codes[$i] == $all_codes[0] ) {
			$same = ($same + 1);
		    }
		}

		if ( $output_selection eq 'VERBOSE' ) {
		    print("same = $same, compare to $Nstr\n");
		    print("neg count = $neg_count\n");
		}

		if ( $same == $Nstr ) {
    
		    if ( $output_selection eq 'VERBOSE' ) {
			print("Not polymorphic in the CC: skip this marker\n\n");
		    }

		} else {

#  include this marker in data file

		    $NCC = ( $NCC + 1 );

		    # ignore chr number in $row[0] 
		    my $bppos = $row[1];
		    my $ref_allele = $row[2];

#  make allele SDP

		    my (@alleles, @transl);

		    $transl[0] = $ref_allele;
		    for ($i=0; $i<$Nalt; $i++) {
			$transl[$i+1] = $alt[$i];
		    }

		    if ( $output_selection eq 'VERBOSE' ) {
                        print("bp position = $row[1]\n");
                        print("allele coding vector @transl\n");
                    }

		    if ($marker_selection eq 'CALL') {

			for ($i=0; $i<=$#all_codes; $i++) {

			    if ( $qual_codes[$i] == -7 ) {
                                my $alc1 = substr( $all_codes[$i], 0, 1);
                                my $alc2 = substr( $all_codes[$i], 2, 1);
                                my $het_str = join('', $transl[$alc1], '/', $transl[$alc2]);
                                $alleles[$i] = $het_str;
#                                print(" HET  $all_codes[$i] qual code $qual_codes[$i] alleles $alc1, $alc2   het string = $het_str\n");
			    } elsif ( $all_codes[$i] >= 0 ) {
				$alleles[$i] = $transl[ $all_codes[$i] ];
			    } elsif ( $all_codes[$i] == -8 ) {
				$alleles[$i] = $header[ $CC_indices[$i] ];
			    } else {
				print("ERROR allele code $all_codes[$i] not recognised\n");
				exit;
			    }
			}

			if ( $output_selection eq 'VERBOSE' ) {
			    print("allele codes @all_codes\n");
			    print("alleles      @alleles\n\n");
			}

		    } elsif ($marker_selection eq 'PRIV') {
			
			for ($i=0; $i<=$#all_codes; $i++) {
			    if ( $all_codes[$i] >= 0 ) {
				$alleles[$i] = $transl[ $all_codes[$i] ];
			    } elsif ( $all_codes[$i] == -8 ) {
				$alleles[$i] = $header[ $CC_indices[$i] ];
			    } else {
				print("ERROR allele code $all_codes[$i] not recognised\n");
				exit;
			    }
			}

			if ( $output_selection eq 'VERBOSE' ) {
			    print("allele codes @all_codes\n");
			    print("alleles      @alleles\n\n");
			}

		    } else {
			print("marker selection $marker_selection not recognised: cannot create allele SDP vector\n");
			exit;
		    }

		    my $out_line = join("\t", $chr, $bppos, @alleles);
		    print( OUTF "$out_line\n");

		} # end of skip marker if loop
	    } # end of skip marker if loop
	} # end of skip marker if loop

#   my @matches = grep { $_ == $bppos } keys %QTL;
#   $QTL{ $bppos }{ 'GeneName' } = [ $gene_name ];
#   $QTL{ $bppos }{ $header[$index] } = [ $ref_allele, 'NA' ];
#    if ( exists ${$QTL{$snppos}}{'Annotation'} ) { 

    } # end of while loop

    close(INF1);
    close(OUTF);

    print("Total number of markers = $Ntotal, number of markers included for CC = $NCC\n");
}
