#! /usr/local/bin/perl

#  map_recomb.perl: calculate number of recombinations per linkage
#  group and individual.  Takes haplotype file as argument in mapmaker
#  format, more or less

# Steve DiFazio, 4-2003

if(!$ARGV[0] ) { die "
_______________________________________________________________\n
Usage: perl map_recomb.perl -a <inputfile> > out.txt \n where
\t-a results in bootstraps including all loci, not just 1 locus per bin,
\t<inputfile> is input file with all genotype info and 5 columns of other data after locus,lg, and pos.
\t<out.txt> is the name of the file to save the output\n
_______________________________________________________________\n\n"; }
warn "\n";

use Getopt::Long;

GetOptions( "a" => \$allloc);

while(<>) {
    chomp();
    if(/^\#/) {
#	($lg,$loc,@prog)=split(/\t/);
	($a,$b,$c,$d,$e,$f,$g,$h,@prog)=split(/\t/);
	
    }
    else {
	my ($loc,$lg,$pos,$obvis,$obinvis,$expectvis,$chi,$diff,@genos)=split(/\t/);
#	($lg,$loc,@genos) = split(/\t/);
	if($#genos ne $#prog) { die "ERROR: $#genos genos and $#prog progeny at line $. of $ARGV\n"; } 
	foreach $i (0 .. $#genos) {
	    if($genos[$i] eq "-" || $counted{$i.":".$lg.":".$pos}) {
		next;
	    }
	    else {
		$genotypes{$genos[$i]} = 1;
		if(!($allloc)) {	$counted{$i.":".$lg.":".$pos}=1; }
		$markers{$lg}[$i]++;
		$totmarkers[$i]++;
		if($genos[$i] eq "H") { 
		    $tricho{$lg}[$i]++; 
		    $tottricho[$i]++;
		}
		if($prevgeno{$lg}[$i]) {
		    if($genos[$i] ne $prevgeno{$lg}[$i]) {
			$block{$genos[$i].":".$lg}[$i] += ($pos - $prevpos{$lg}[$i])/2;
			$block{$prevgeno[$i].":".$lg}[$i] += ($pos - $prevpos{$lg}[$i])/2;
			$block{$genos[$i]}[$i] += ($pos - $prevpos{$lg}[$i])/2;
			$block{$prevgeno[$i]}[$i] += ($pos - $prevpos{$lg}[$i])/2;
			$recomb{$lg}[$i]++;
			$totrecomb[$i]++;
			$prevgeno{$lg}[$i] = $genos[$i];
		    }
		    else {
			$cblock{$genos[$i].":".$lg}[$i] += ($pos - $prevpos{$lg}[$i]);
			$cblock{$genos[$i]}[$i] += ($pos - $prevpos{$lg}[$i]);
			$block{$genos[$i].":".$lg}[$i] += ($pos - $prevpos{$lg}[$i]);
			$block{$genos[$i]}[$i] += ($pos - $prevpos{$lg}[$i]);
		    }
		}
		else{
		    $prevgeno{$lg}[$i] = $genos[$i];
		    $block{$genos[$i].":".$lg}[$i] += ($pos - $prevpos{$lg}[$i]);
		    $block{$genos[$i]}[$i] += ($pos - $prevpos{$lg}[$i]);
		}
		$prevpos{$lg}[$i] = $pos;
	    }
	}
    }
}
print "LG\tData\t";
foreach $i (0 .. $#prog) {
    print "$prog[$i]\t";
}
print "\n";
foreach $lg (sort (keys(%prevgeno))) {
    print "$lg\tMarkers\t";
    foreach $i (0 .. $#prog) {
	print $markers{$lg}[$i] ? $markers{$lg}[$i] : 0,"\t";
    }
    print "\n";
    print "$lg\tRecombinations\t";
    foreach $i (0 .. $#prog) {
	print $recomb{$lg}[$i] ? $recomb{$lg}[$i] : 0,"\t";
    }
    print "\n";
    print "$lg\tTricho\t";
    foreach $i (0 .. $#prog) {
	print $tricho{$lg}[$i] ? $tricho{$lg}[$i] : 0,"\t";
    }
    print "\n";
    
    foreach $geno (keys(%genotypes)) {
	print "$lg\t$geno block\t";
	foreach $i (0 .. $#prog) {
	    print $block{$geno.":".$lg}[$i] ? $block{$geno.":".$lg}[$i] : 0,"\t";
	}
	print "\n";
	print "$lg\t$geno cblock\t";
	foreach $i (0 .. $#prog) {
	    print $cblock{$geno.":".$lg}[$i] ? $cblock{$geno.":".$lg}[$i] : 0,"\t";
	}
	print "\n";
    }
}
    

print "Total\tMarkers\t";
foreach $i (0 .. $#prog) {
    print $totmarkers[$i] ? $totmarkers[$i] : 0,"\t";
}
print "\n";
print "Total\tRecombinations\t";
foreach $i (0 .. $#prog) {
    print $totrecomb[$i] ? $totrecomb[$i] : 0,"\t";
}
print "\n";
print "Total\tTricho\t";
foreach $i (0 .. $#prog) {
    print $tottricho[$i] ? $tottricho[$i] : 0,"\t";
}
print "\n";
foreach $geno (keys(%genotypes)) {
    print "Total\t$geno block\t";
    foreach $i (0 .. $#prog) {
	print $block{$geno}[$i] ? $block{$geno}[$i] : 0,"\t";
    }
    print "\n";
    print "Total\t$geno cblock\t";
    foreach $i (0 .. $#prog) {
	print $cblock{$geno}[$i] ? $cblock{$geno}[$i] : 0,"\t";
    }
    print "\n";
}

sub by_num {
    $a <=> $b;
}