#!/usr/bin/perl

use strict;

my $dir1 = "/Net/mus/data/www/preCC/MUGA/GENOTYPES/";
my $dir2 = "/Net/mus/data/www/preCC/MEGA_MUGA/Feb2014.MEGA+MDA/";
my $dir3 = "/Net/mus/data/www/preCC/MEGA_MUGA/Feb2015/gstemp/";

my $outdir = "/Net/mus/data/www/preCC/MEGA_MUGA/Mar2015.MEGA+MDA+MUGA";
my $prefix = "MEGA+MDA+MUGA";
my $fam1 = "MUGA";
my $fam2 = "MEGA+MDA";

mkdir( $outdir );

my @strain1 = ("AJ", "C57BL6J", "129S1SvImJ", "NODLtJ", "NZOHlLtJ", "CASTEiJ", "PWKPhJ", "WSBEiJ");
my @strain2 = ( "A/J", "C57BL/6J", "129S1SvImJ", "NOD/LtJ", "NZO/HILtJ", "CAST/EiJ", "PWK/PhJ", "WSB/EiJ");
my $strain;
for ( my $k=0;$k<8;$k++) {
  $strain->{$strain1[$k]} = $strain2[$k];
}

my $istrain;
foreach my $s ( keys %$strain ) {
  $istrain->{$strain->{$s}} = $s;
}

my $csv = readCSVGenotypes( $dir3);
my $map3 = $csv->{map};
my $nchr = 20;


open(DEBUG, ">debug.txt");

for ( my $c=1;$c<=$nchr;$c++) {
  my  $intersect;
  my $m=0;
  my @map1;
  my @map2;
  my $map1;
  my $map2;
  my $snp1;
  my $snp2;
  my $allele1;
  my $allele2;
  my $include;
  my $exclude;

  my $chr = ( $c == 20 ) ? "X" : $c;

  my $fil1 = "$dir1/chr$chr.muga.map";
  open( MAP, $fil1 ) ||die "Could not open $fil1\n";
  $_ = <MAP>;
  while (<MAP>) {
    chomp;
    my ($snp, $ch, $bp ) = split;
    $ch =~ s/chr//;
    $map1->{$chr}{$bp} = $snp;
    push @map1, [$snp, $chr, $bp ];
    $snp1->{$snp} = [$chr, $bp ];
  }
  print DEBUG "MUGA $fil1 $#map1\n";
  
  my $fil2 = "$dir2/chr$chr.MEGA+MDA.map";
  open( MAP, $fil2 ) ||die "Could not open $fil2\n";
  $_ = <MAP>;
  while (<MAP>) {
    chomp;
    my ($snp, $ch, $bp ) = split;
    $ch =~ s/chr//;
    $map2->{$chr}{$bp} = $snp;
    push @map2, [$snp, $chr, $bp ];
    $snp2->{$snp} = [$chr, $bp ];
  }
  print DEBUG "MEGA+MDA $fil2 $#map2\n";

  
  my $n=0;
  foreach my $bp1 ( sort {$a<=>$b} keys %{$map1->{$chr}} ) {
    #    print "$chr $bp1\n";
    if ( ($map2->{$chr}{$bp1} eq $map2->{$chr}{$bp1}) && exists( $map3->{$chr}{$bp1} ) ) {
      $intersect->{$chr}{$bp1} = 1;
      $n++;
    }
  }

  print DEBUG "intersect $chr $n $m\n";

  $m += $n;
  print "$chr $n $m\n";

 
  my $fil1 = "$dir1/chr$chr.muga.alleles";
  open( ALLELES, $fil1 ) ||die "Could not open $fil1\n";
  $_ = <ALLELES>;
  chomp;
  
  my ( $x, $nm, $y, $ns ) = split;
  $_ = <ALLELES>;
  chomp;
  my ( $x, @s ) = split;
  my $snp;
  while (<ALLELES>) {
    if ( /^marker/ ) {
      my @a = split;
      $snp = $a[1];
    } else {
      chomp;
      my ( $x, $a, @a ) = split;
      my $k=0;
      if ( $a ne "NA" ) {
	foreach my $s ( @strain1 ) {
	  $allele1->{$snp}{$strain->{$s}} = $a if ( $a[$k]>0 );
	  $k++;
	}
      }
    }
  } 
  close(ALLELES);

  my $fil2 = "$dir2/chr$chr.MEGA+MDA.alleles";
  open( ALLELES, $fil2 ) ||die "Could not open $fil2\n";
  $_ = <ALLELES>;
  chomp;
  
  my ( $x, $nm, $y, $ns ) = split;
  $_ = <ALLELES>;
  my $snp;
  while (<ALLELES>) {
    if ( /^marker/ ) {
      my @a = split;
      $snp = $a[1];
    } else {
      chomp;
      my ( $x, $a, @a ) = split;
      my $k=0;
      if ( $a ne "NA" ) {
	foreach my $s ( @strain2 ) {
	  $allele2->{$snp}{$s} = $a if ( $a[$k]>0 );
	  $k++;
	}
      }
    }
  }
  close(ALLELES);

  my $consistent=0;
  my $inconsistent=0;
  foreach my $bp ( keys %{$map1->{$chr}} ) {
    my $snp1 = $map1->{$chr}{$bp};
    if ( exists( $map2->{$chr}{$bp} ) ) {
      my $snp2 = $map2->{$chr}{$bp};
      my $y=0;
      foreach my $s ( keys %{$allele1->{$snp1}} ) {
	$y += $allele1->{$snp1}{$s} eq $allele2->{$snp2}{$s};
      }
      print "$chr $bp $snp1 $snp2 $y\n";
      if ( $y==8 ) {
	$include->{$snp1}{$snp2} = ( $y==8 );
	$consistent++;
      }
      else {
	$inconsistent++;
      }
    }
  }
  print DEBUG "$chr consistent $consistent inconsistent $inconsistent\n";
  my $excluded=0;
  my $included=1;
  foreach my $snp1 ( keys %{$map1->{$chr}} ) {
    if ( ( exists( $map2->{$chr}{$snp1}) && ! exists( $include->{$snp1}{$snp1} ) ) ) {
      $exclude->{$snp1} = 1;
      $excluded++;
    }
    else {
      $included++;
    }
  }

  print DEBUG "$chr excluded $excluded included $included\n";
  
  my $n1=0;
  my $n2=0;
  my $bp1;
  my $bp2;
  my @combined;
  print "here 2 $#map1 $#map2\n";
  $exclude->{"JAX00712235"} = 1;
  
  do {
    my $ref1 = $map1[$n1];
    my $ref2 = $map2[$n2];
    my $snp1 = $ref1->[0];
    my $snp2 = $ref2->[0];
    $bp1 = $ref1->[2];
    $bp2 = $ref2->[2];
    #    print "$n1 $n2 $bp1 $bp2\n";
    if ( $bp1 < $bp2 ) {
      if ( ! $exclude->{$snp1} ) {
	push @combined, [$snp1, $chr, $bp1, 1, undef];
      }
      $n1++;
    } elsif ( $bp1 == $bp2 ) {
      if ( $include->{$snp1}{$snp2} ) {
	push @combined, [$snp1, $chr, $bp1, 3, $snp2 ];
      }
      $n1++;
      $n2++;
    } elsif ( $bp1 > $bp2 ) {
      if ( ! $exclude->{$snp2}) {
	push @combined, [$snp2, $chr, $bp2, 2, undef];
      }
      $n2++;
    }
  } while ( $n1 <= $#map1 && $n2 <= $#map2 && $bp1 && $bp2 );
  my $nm = $#combined+1;
  print DEBUG "$chr combined $nm\n";
  
  print "here 2\n";

  my $mapfile = "$outdir/chr$chr.$prefix.map";
  open(MAP, ">$mapfile") || die "Could not write $mapfile\n";
  print MAP "marker	chromosome	bp\n";
  foreach my $ref ( @combined ) {
    print MAP join( "\t", $ref->[0], $ref->[1], $ref->[2] ) . "\n";
  }
  close(MAP);
  print "here 1\n";
  my $allelesfile = "$outdir/chr$chr.$prefix.alleles";
  open(ALLELES, ">$allelesfile") || die "Could not write $allelesfile\n";
  print ALLELES "markers $nm strains 8\nstrain_names @strain2\n";
  my $allele;
  my $hetg = {};
  foreach my $ref ( @combined ) {
    if ( $ref->[3] == 1 || $ref->[3] == 3 ) {
      $allele = $allele1;
    } else {

      $allele = $allele2;
    }
    my $cm = $ref->[2]* 5 * 1e-7;
    my $snp = $ref->[0];
    print ALLELES "marker $snp 3 $ref->[1] $cm $ref->[3]\n";
    my %is_a;
    foreach my $s ( keys %{$allele->{$snp}} ) {
      $is_a{$allele->{$snp}{$s}}++;
    }
    foreach my $a ( sort keys  %is_a ) {
      my $x = 1/$is_a{$a};
      print ALLELES "allele $a";
      foreach my $s ( @strain2 ) {
	printf ALLELES "\t%.4f", $allele->{$snp}{$s} eq $a ? $x: 0;
      }
      print ALLELES "\n";
    }
    $hetg->{$ref->[1]}{$ref->[2]} = join(  " ", sort keys  %is_a );
#    print DEBUG "hetg $ref->[1] $ref->[2]  $hetg->{$ref->[1]}{$ref->[2]}\n";
    print ALLELES "allele NA";
    for (my $k=1;$k<=8;$k++) {
      print ALLELES "\t0.1250";
    }
    print ALLELES "\n";
    
  }
  
  close(ALLELES);
  
  my $ped1 = "$dir1/chr$chr.muga.data"; 
  my $ped2 = "$dir2/chr$chr.MEGA+MDA.data";
  my $ped  = "$outdir/chr$chr.$prefix.data";
  
  open(PED1,$ped1) || die "Could not open $ped1\n";
  open(PED2,$ped2) || die "Could not open $ped2\n";
  open(PED,">$ped") || die "Could not open $ped\n";
  my $count;
  while (<PED1>) {
    chomp;
    my $g;
    my ( $fam, $id, $dad, $mum, $sex, $p, @g ) = split;
    $id =~ s/IL\./IL-/;
    for ( my $k=0, $m=0;$k<$#g;$k+=2,$m++) {
      $g->{$map1[$m]->[0]} = $g[$k] . " " . $g[$k+1];
    }
    $count->{$id}++;
    print PED "$fam1 $id $dad $mum $sex $p";
    foreach my $ref (@combined ) {
      my $snp = $ref->[0];
      my $x = defined($g->{$snp}) ? $g->{$snp} : "NA NA";
      print PED "\t$x";  
    }
    print PED "\n";
  }
  close(PED1);
  while (<PED2>) {
    chomp;
    my $g;
    my ( $fam, $id, $dad, $mum, $sex, $p, @g ) = split;
    $id =~ s/IL\./IL-/;
    if ( $count->{$id} ) {
      $fam = "$fam";
      $id = "$id-$fam";
    }
    $count->{$id}++;
    for ( my $k=0, $m=0;$k<$#g;$k+=2,$m++) {
      $g->{$map2[$m]->[0]} = $g[$k] . " " . $g[$k+1];
    }
    
    print PED "$fam $id $dad $mum $sex $p";
    
    foreach my $ref (@combined ) {
      my $snp = ( $ref->[3] == 3 ) ? $ref->[4] : $ref->[0];
      my $x = defined($g->{$snp}) ? $g->{$snp} : "NA NA";
      print PED "\t$x";  
    }
    print PED "\n";
  }

  close(PED2);
  my $csvg = $csv->{data};
  my $k=0;
  foreach my $id ( sort keys %$csvg )  {
    my $csvgid = $csvg->{$id};
    $id = "$id-MEGA.UNC" if ( $count->{$id} );
      
    print PED "MEGA.UNC $id 0 0 1 0";
    foreach my $ref (@combined ) {
      my $g = $csvgid->{$ref->[1]}{$ref->[2]};
      if ( ($g eq "N") || (!$g) ) {
	$g = "NA NA";
      } elsif ( $g eq "H" ) {
	$g = $hetg->{$ref->[1]}{$ref->[2]};
      }
      else { $g = "$g $g" }
      print PED " $g";
     print DEBUG  "$ref->[1]|$ref->[2]|$csvgid->{$ref->[1]}{$ref->[2]}|$g|\n" if ( !$g ) ;
    }
    print PED "\n";
  }
	
    close(PED);
}



sub readCSVGenotypes {
    my ( $dir, $max ) = @_;

    my @files = glob( "$dir/*.csv" );
    my $data = {};
    my $map = {};
    @files = splice( @files, 0, $max ) if ( $max > 0 ) ;
    foreach my $fil ( @files ) {
      my ( $id ) = $fil =~ /\/(\w+)_\w+.csv$/;
      $id =~ s/IL/IL-/;
      open(FILE, $fil) || die "Coudl not open $fil\n";
      $_ = <FILE>;
      my $n=0;
      while (<FILE>) {
	chomp;
	my ($chr,$bp,$snp,$f) = split(/\,/);
	$map->{$chr}{$bp}++;
	$data->{$id}{$chr}{$bp} = $snp;
	$n++;
      }
      close(FILE);
      print "$fil $id $n\n";
    }

    return ( {map=>$map, data=>$data} );
  }
