
#!/usr/bin/perl

use strict;


my $fname = './final.vcf';

my $file_chr = $ARGV[0];
my $outfile = join('', 'final.chr', $file_chr, '.vcf.txt');


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

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

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

while ( $start eq '##' ) {

    $next = substr( $line, 2, 10);

    if ( $next eq 'FORMAT=ATG' ) {
# write out string to outfile
	print( OUTF "$line");
    }

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

my $header_line = substr( $line, 1);
print("$header_line\n");
my @header = split(/\t/, $header_line);

my @indices = (0, 1, 3, 4, 9..25);
my @new_header = @header[@indices];
my $header_string = join("\t", @new_header);
print( OUTF "$header_string");


my ($line, $chr, $bppos, $ref, $alt, $str);
my (@row, @new_row, @info_vec, @strains, @strain_info, @geno_vec);

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

    @row = split(/\t/, $line); 
    @new_row = @row[@indices];

    $chr = $new_row[0];

    if ( $file_chr == $chr ) {

	$bppos = $new_row[1];
	$ref = $new_row[2];
	$alt = $new_row[3];

	@strains = @new_row[ 4..$#new_row ];

	my @geno_codes;
	my @geno_quals;
	my @out_data;

	foreach $str ( @strains ) {

	    @strain_info = split(/:/, $str);

	    my $geno;
	    @geno_vec = split(/\//, $strain_info[0]);

	    if ( $#geno_vec == 1 ) {
		if ( $geno_vec[0] == $geno_vec[1] ) {
		    $geno = $geno_vec[0];
		} else {
		    if ( $strain_info[1] == -7 ) {
			$geno = $strain_info[0];
		    } else {
			print("what is this genotype?  @geno_vec   @strain_info\n");
			exit;
		    }
		}
	    } else {
		print("ERROR unexpected genotype string @geno_vec qual $strain_info[1]\n");
		exit;
	    }

	    push( @geno_codes, $geno );
	    push( @geno_quals, $strain_info[1] );

	    my $data = join(':', $geno, $strain_info[1] );
	    push( @out_data, $data );
	}

#	print("genotype quality scores: @geno_quals\n");
#	print("genotype codes: @geno_codes\n");
#	print("output = @out_data\n\n");

	my $out_vec = join("\t", $chr, $bppos, $ref, $alt, @out_data);
	print( OUTF "$out_vec\n");

    } # end of chr 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(OUTF);

close(INF1);
