#!/usr/bin/perl

use strict;


my $output_selection = 'CONCISE'; # CONCISE or VERBOSE
my $marker_selection = 'CALL'; # 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_HS = 1..19;
push (@chrs_HS, 'X');

my @HS_indices = (8, 7, 9, 10, 11, 13, 14, 15);

foreach my $chr ( @chrs_HS ) {

    my $fname = join('', 'rawdata.alleles.17strains.chr', $chr);
    print("\n$fname\n");

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


    my $i;
    my $Ntotal = 0;
    my $NHS = 0;

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

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

    my @HS_strains = @header[@HS_indices];
    my $out_header = join("\t", 'chr', 'bp', @HS_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 HS strains

	my @all_codes;
	foreach my $index ( @HS_indices ) {
	    push(@all_codes, $row[$index]);
	}

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

#  is the marker polymorphic in the HS strains?

	my $pos_count = 0;
	my $neg_count = 0;
	my $unc_calls = 0;
	my $same = 0;
	my $Nstr = scalar @all_codes;

	for ($i=0; $i<=$#all_codes; $i++) {
	    if ( $all_codes[$i] > 0 ) {
		$pos_count = ($pos_count + 1);
	    } elsif ( $all_codes[$i] < 0 ) {
		$neg_count = ($neg_count + 1);
		if ( $all_codes[$i] < -1 && $all_codes[$i] > -8 ) {
		    $unc_calls = ($unc_calls + 1);
		}
	    }
	    if ( $all_codes[$i] == $all_codes[0] ) {
		$same = ($same + 1);
	    }
	}

	if ( $output_selection eq 'VERBOSE' ) {
	    print("same = $same, compare to $Nstr\n");
	    print("pos count = $pos_count\n");
	    print("neg count = $neg_count\n");
	    print("uncertain calls for multi-allelic markers = $unc_calls\n");
	}

	if ( $same == $Nstr ) {

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

        } else {

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

            my $alt_allele = $row[3];
            my $Nalt = length( $alt_allele );

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

	    if ( $marker_selection == '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 values

		if ($marker_selection eq 'CALL') {

		    if ( $neg_count > 0 ) {
			for ($i=0; $i<= $#all_codes; $i++) {
			    if ( $all_codes[$i] == -1 ) {
				# call as ref
				$all_codes[$i] = 0;
			    } elsif ( $all_codes[$i] == -8 ) {
				# leave as neg score
			    } elsif ( $all_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
				$all_codes[$i] = 1;
			    }
			}
		    } # end of neg count if loop

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

#  is the marker still polymorphic in the HS strains?

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

		for ($i=0; $i<=$#all_codes; $i++) {
		    if ( $all_codes[$i] > 0 ) {
			$pos_count = ($pos_count + 1);
		    } elsif ( $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("pos count = $pos_count\n");
		    print("neg count = $neg_count\n");
		}

		if ( $same == $Nstr ) {

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

		} else {

#  include this marker in data file

		    $NHS = ( $NHS + 1 );

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

#  make allele coding vector

		    my @all;
		    $all[0] = $ref_allele;
                
		    if ( $Nalt > 0 ) {
			for ($i=0; $i<$Nalt; $i++) {
			    $all[($i+1)] = substr( $alt_allele, $i, 1);
			}
		    }

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

#  make allele SDP

		    my @alleles;

		    if ($marker_selection eq 'CALL') {

			my $het_str = join('', $all[0], '/', $all[1]);

			for ($i=0; $i<=$#all_codes; $i++) {
			    if ( $all_codes[$i] >= 0 ) {
				$alleles[$i] = $all[ $all_codes[$i] ];
			    } elsif ( $all_codes[$i] == -8 ) {
				$alleles[$i] = $header[ $HS_indices[$i] ];
			    } elsif ( $all_codes[$i] == -7 ) {
				$alleles[$i] = $het_str;
			    } else {
				print("ERROR allele code $all_codes[$i] not recognised\n");
				exit;
			    }
			}

			if ( $output_selection eq 'VERBOSE' ) {
			    print("het string $het_str\n");
			    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] = $all[ $all_codes[$i] ];
			    } elsif ( $all_codes[$i] < 0 && $all_codes[$i] >= -8 ) {
				$alleles[$i] = $header[ $HS_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 HS = $NHS\n");
}
