#!/usr/bin/perl

use strict;

my $build37 = "build37.map";
my $mm10 = "hglft_genome_751b_dcc510.bed";

open(BUILD37,$build37) || die "Could not open $build37\n";

my @b37;
my $snp37;
while(<BUILD37>) {
    chomp;
    my @a = split;
    push @b37, \@a;
    $snp37->{$a[3]} = \@a;
}
close(BUILD37);

open(MM10,$mm10) || die "Could not open $mm10\n";

my @m10;
my $snp10;
while(<MM10>) {
    chomp;
    my @a = split;
    push @m10, \@a;
    $snp10->{$a[3]} = \@a;
}
close(MM10);

my @map = sort  { my $x = $a->[0] cmp $b->[0]; if ( !$x ) { $x = $a->[1]-$b->[1] } return($x) } @m10;
my $chrmap = {};
my $chrsnp = {};
my $chrorder = {};
my $nsnp = {};
foreach my $ref (@map) {
    my $chr = $ref->[0];
    $chrmap->{$chr} = () if ( !exists($chrmap->{$chr}));
    $chrorder->{$chr} = () if ( !exists($chrorder->{$chr}));
    $nsnp->{$chr} = 0 if ( !exists($nsnp->{$chr}));
    push @{$chrmap->{$chr}}, $ref;
    $chrsnp->{$ref->[3]} = $nsnp->{$chr};
    $chrorder->{$chr}[$nsnp->{$chr}] = $ref->[3];
    $nsnp->{$chr}++;
}

for( my $chr=1;$chr<=20;$chr++) {
    my $c =  ( $chr < 20) ? "chr$chr" : "chrX"; 
    my $afile = "$c.MEGA+MDA+MUGA.alleles";
    if ( ! -s $afile ) {
	my $cmd = "curl http://mtweb.cs.ucl.ac.uk/mus/www/preCC/MEGA_MUGA/Mar2015.MEGA+MDA+MUGA/Mar2015.MEGA+MDA+MUGA/chr$c.MEGA+MDA+MUGA.alleles > $afile";
	system($cmd);
	print "$cmd\n";
    }
    
    open( AFILE, $afile ) || die "Could not open $afile\n";
    $_ = <AFILE>;
    my $snames = <AFILE>;
    my $atxt = {};
    my ( $id, $na, $ch, $cm, $bp );
    my $m=0;
    my $ipos = {};
    while(<AFILE>) {
	if ( /^marker\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)/ ) {
	    ( $id, $na, $ch, $cm, $bp ) = ($1, $2, $3, $4, $5);
	    my $ref = $snp10->{$id};
	    if ( $ref) {
		my $bp10 = $ref->[2];
		my $cm10 = 1500*$bp10/3.0e9;
		$atxt->{$id} = "marker $id $na $ch $cm10 $bp10\n";
	    }
	    else {
		$atxt->{$id} = "$_\n";
	    }
	    $ipos->{$id} = $m++;
	} else {
	    $atxt->{$id} .= $_;
	}
    }
    close(AFILE);

    my $afile = "./MM10/$afile\n";
    open(AFILE, ">$afile") || die "Could not write $afile\n";
    print AFILE "markers $nsnp->{$c} strains 8\n$snames";
    for( my $k=0;$k<$nsnp->{$c};$k++) {
	print AFILE $atxt->{$chrorder->{$c}[$k]};
	die "no data for snp $chr $k $chrorder->{$c}[$k]\n" if ( ! $atxt->{$chrorder->{$c}[$k]} );
    }
    close(AFILE);

    my $pedfile = "$c.MEGA+MDA+MUGA.data";

    if ( ! -s $pedfile ) {
	my $cmd = "curl http://mtweb.cs.ucl.ac.uk/mus/www/preCC/MEGA_MUGA/Mar2015.MEGA+MDA+MUGA/Mar2015.MEGA+MDA+MUGA/chr$c.MEGA+MDA+MUGA.data > chr$c.MEGA+MDA+MUGA.data";
	system($cmd);
	print "$cmd\n";
    }

    open(PED, $pedfile) || die "Could not open $pedfile\n";
    my $pedout = "MM10/$pedfile";
    open(PEDOUT, ">$pedout") || die "Could not write $pedout\n";

    my @data;
    while(<PED>) {
	chomp;
	my ( $fam, $id, $dad, $mum, $sex, $phen, @a ) = split;
	my @g;
	for( my $k=0, my $m=0;$k<$#a;$k+=2,$m++) {
	    $g[$m] =  $a[$k] . " " . $a[$k+1];
	}
	my @h;
	for( my $k=0;$k<$nsnp->{$c};$k++) {
	    push @h, $g[$ipos->{$chrorder->{$c}[$k]}];
	}
	print PEDOUT join( "\t", $fam, $id, $dad, $mum, $sex, $phen, @h) . "\n";
    }
    close(PED);
    close(PEDOUT);

}

	
