#!/usr/bin/perl

use strict;

my @files = @ARGV;
my @founders = ( "A/J","C57Bl/6J","129SvlmJ","NOD/LtJ","NZO/HlLtJ","CAST/EiJ","PWK/PhJ","WSB/EiJ");

my $ok = {"A"=>" A A", "C"=>" C C", "G"=>" G G", "T"=>" T T", "K"=>" A T", "M"=>" A C", "Y"=>" C T", "R"=>" A G", "H"=>"H"};
my $het;
my $error = 0.001;

foreach my $f (@files) {
  print "reading $f\n";
  my $ccline;
  my $ccidx;
  my $snp;
  my $hallele;
    
  my ($chr) = $f =~ /chr(.+).csv$/;
  $chr = $chr+0 if ( $chr =~ /^\d+$/);
	       
  my $pedfile = "chr$chr.preCC.ped";
  my $allelesfile = "chr$chr.preCC.alleles";
  my $mapfile = "chr$chr.preCC.map";

  open(FILE,$f) || die "could not open $f\n";
  $_ = <FILE>;
  chomp;
  my @header = split(/,/);
  my $colidx;
  for (my $k=0;$k<=$#header;$k++) {
    $colidx->{$header[$k]} = $k;
    if ( $header[$k] =~ /^(OR\d+)/ ) {
      $ccline->{$header[$k]} = $header[$k];
      $ccidx->{$header[$k]} = $k;
    }
  }

  my $n=0;
  my ( $geno, $fgeno, $fallele, $hgeno );
  while (<FILE>) {
    chomp;
    my @a = split(/,/);
    $snp->[$n] = [$a[$colidx->{id}], $a[$colidx->{chr}], $a[$colidx->{position}], $a[$colidx->{call0}], $a[$colidx->{call2}], " $a[$colidx->{call0}] $a[$colidx->{call2}]"  ];
    $hallele->[$n] = " $a[$colidx->{call0}] $a[$colidx->{call2}]";
    $ok->{H} = " $a[$colidx->{call0}] $a[$colidx->{call2}]";
    for (my $k=6;$k<14;$k++) {
      $fgeno->[$n]{$header[$k]} = $a[$k];
      #	    print "$n $k,$a[$k],$ok->{$a[$k]}\n" if ( $n < 10 );
    }
    for (my $k=14;$k<=$#a;$k++) {
      $geno->[$n]{$header[$k]} = $ok->{$a[$k]};
      $hgeno->[$n]{$header[$k]} = $a[$k];

    }

    my $allele;
    foreach my $p (@founders ) {
      $allele->{$fgeno->[$n]{$p}}++;
    }
    my $na=0;
    my $nh=0;
    my $miss=0;
    #	print join( ",", ">> $n", sort keys %$allele) . "\n" if ( $n < 10 );
    foreach my $a ( sort keys %$allele ) {

      if ( $a eq "H" ) {
	$nh++;
      } elsif ( $a =~ /^[\sACGT]$/ ) {
	$na++;
      } else {
	#	    print "$n alleles $a\n" if ( $n < 100 ) ;
	$miss++;
      }
    }

    $fallele->[$n] = {allele=>$allele, na=>$na, miss=>$miss, het=>$nh};
    if ( $n < 100 ) {
      my @k = keys %$allele;
      #	  print "$n $allele @k $na $miss $nh\n";
    }

    $n++;
  }


  my $ns = $#$snp +1;
  my $good=0;
  my $use;

  for (my $n=0;$n<$ns;$n++) {
    #      print "$n $good $fallele->[$n]{miss} $fallele->[$n]{na} $fallele->[$n]{het}\n";
    if ( $fallele->[$n]{miss} == 0 && $fallele->[$n]{na} == 2 && $fallele->[$n]{het} ==0 ) {
      $use->[$n] = 1;
      $good++;
      my @v = sort keys %{$fallele->[$n]{allele}};

      #	print "$n $hallele->[$n]\n" ;
    }
  }
    
  open(PED, ">$pedfile") || die "could not write $pedfile\n";
  print "writing $pedfile\n";
  foreach my $line ( sort keys %$ccline ) {
    printf PED "preCC\t%s\t0\t0\tM\tNA", $ccline->{$line};
    for (my $k=0; $k<=$#$snp;$k++) {
      if ( $use->[$k] ) {
	my $g = $geno->[$k]{$line};
	$g = $g eq "H" ? $hallele->[$k] : $g;
	print PED $g ? $g : " NA NA";
	$het->{$line}{$hgeno->[$k]{$line}}++;
      }
    }
    print PED "\n";
  }
  close(PED);
    
  foreach my $line ( keys %$het ) {
    my $tot=0;
    my $h = 0;
    my $n = 0;
    foreach my $a ( keys %{$het->{$line}} ) {
      print "$a $het->{$line}{$a}\n";
      $tot += $het->{$line}{$a};
      $n += $het->{$line}{$a} if ( $het->{$line}{$a} eq "N");
    }
    $h = $het->{$line}{"H"};
    $tot -= $n;
    my $q = $h/$tot;
    print "$line calls $tot $n $h $q\n";
  }

  open(ALLELES, ">$allelesfile\n") || die "could not write $allelesfile\n";
  open(MAP, ">$mapfile\n") || die "could not write $mapfile\n";
  print MAP "marker\tchromosome\tbp\n";
  print ALLELES "markers $good strains 8\nstrain_names " . join( " " , @founders ) . "\n";
  for (my $n=0;$n<$ns;$n++) {
    if ( $use->[$n] ) {
      print MAP join( "\t", $snp->[$n][0], $chr , $snp->[$n][2]) . "\n";
      my $allele = $fallele->[$n]{allele};
      my $na = $fallele->[$n]{na};
      my $na1=$na+1;
      my $cm = $snp->[$n][2]/2e6;
      print ALLELES "marker $snp->[$n][0] $na1 $chr $cm\n";
      foreach my $a ( sort keys %$allele ) {
	if ( $a =~ /[ACGT]/ ) {
	  print ALLELES "allele $a ";
	  foreach my $s ( @founders ) {
	    my $p = $fgeno->[$n]{$s} eq $a ? 1/$allele->{$a} : $error;
	    printf ALLELES " %.4f", $p;
	  }
	  print ALLELES "\n";
	}
      }
      print ALLELES "allele NA ";
      foreach my $s ( @founders ) {
	print ALLELES " 0.125";
      }
      print ALLELES "\n";
    }
  }
  close(ALLELES);
  close(MAP);
  print "nsnps $good/$ns \n";
}
