#!/usr/bin/perl

use strict;

my $dir1 = "./Feb2014/";
my $dir2 = "/data/www/preCC/preCC.merge/";
my $dir3 = "./Feb2014.MEGA+MDA";
my $prefix = "MEGA+MDA";
my $fam1 = "MDA";
my $fam2 = "MEGA";

mkdir( $dir3 );

my @strain1 = ( "129S1/SvImJm", "A/Jm0417", "C57BL/6Jm38420", "CAST/EiJm0538", "NOD/ShiLtJm1214", "NZO/HILtJm0591", "PWK/PhJm1257", "WSB/EiJm1345" );
my @strain2 = ( "A/J", "C57BL/6J", "129S1SvImJ", "NOD/LtJ", "NZO/HILtJ", "CAST/EiJ", "PWK/PhJ", "WSB/EiJ");
my $strain;
$strain->{"129S1/SvImJm"} = "129S1SvImJ";
$strain->{"A/Jm0417"} = "A/J";
$strain->{"C57BL/6Jm38420"} = "C57BL/6J";
$strain->{"CAST/EiJm0538"} = "CAST/EiJ";
$strain->{"NOD/ShiLtJm1214"} = "NOD/LtJ";
$strain->{"NZO/HILtJm0591"} = "NZO/HILtJ";
$strain->{"PWK/PhJm1257"} =  "PWK/PhJ";
$strain->{"WSB/EiJm1345"}  = "WSB/EiJ";

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

for ( my $c=1;$c<=20;$c++) {
my $chr = $c==20 ? "X" : $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 $fil1 = "$dir1/chr$chr.MEGAMUGA.MERGED.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 ];
  }

  my $fil2 = "$dir2/chr$chr.merged.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 ];
  }

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

 
  my $fil1 = "$dir1/chr$chr.MEGAMUGA.MERGED.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.merged.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);

  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";
      $include->{$snp1}{$snp2} = ( $y==8 );
    }
  }

  foreach my $snp1 ( keys %{$map1->{$chr}} ) {
    $exclude->{$snp1} = ( exists( $map2->{$chr}{$snp1}) && ! exists( $include->{$snp1}{$snp1} ) );
  }

  my $n1=0;
  my $n2=0;
  my $bp1;
  my $bp2;
  my @combined;
  print "here 2 $#map1 $#map2\n";

  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 "here 2\n";

  my $mapfile = "$dir3/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 = "$dir3/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;
  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";
    }
    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.MEGAMUGA.MERGED.data"; 
  my $ped2 = "$dir2/chr$chr.merged.data";
  my $ped  = "$dir3/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} ) {
      for( my $k=0, $m=0;$k<$#g;$k+=2,$m++) {
	$g->{$map2[$m]->[0]} = $g[$k] . " " . $g[$k+1];
      }
      
      print PED "$fam2 $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);
  close(PED);
}
  
