#!/usr/bin/perl -w

# ld50calc.pl
# Purpose: 
# Input: Survival Chart, with format first line: Headers for all columns (strain and day),
#        First column containing doses from HIGHEST TO LOWEST with all doses having same dilution factor,
#        All other columns with number of SURVIVING mice
# Output: Tab-delimited table of LD50 values
# Written by Charlie Kim
# September 27, 2000

#########################
# Introductory Comments #
#########################

print "All files can be specified as command line arguments by using flags.\n";
print "    -i [survival numbers tab file]\n";
print "    -o [output filename]\n";
print "Note that the input file must have a header for each column.\n\n";

##############
# Open Files #
##############


use Getopt::Std;
getopts('i:o:');

if ($opt_i) { open(LDFILE, $opt_i) or die "$opt_i not found!\n"; }
else {
    print 'Enter your tab-delimited LD50 survival data filename: ';
    chomp ($ldfile = <STDIN>);
    open(LDFILE, $ldfile) or die "$ldfile not found!\n";
}

if ($opt_o) {
    if (-e $opt_o) {
	print "$opt_o already exists!  Overwrite (Y/N)? ";
	chomp ($_ = <STDIN>);
	while (/[^yn]/i) {
	    print 'Y or N, please: ';
	    chomp ($_ = <STDIN>);
	}
	if (/n/i) { die "$opt_o not overwritten.\n"; }
	else { open(OUTFILE, ">$opt_o"); }
    }
    else { open(OUTFILE, ">$opt_o"); }
}

else {
    print 'Enter a LD50 output filename: ';
    chomp ($outfile = <STDIN>);
    if (-e $outfile) {
	print "$outfile already exists!  Overwrite (Y/N)? ";
	chomp ($_ = <STDIN>);
	while (/[^yn]/i) {
	    print 'Y or N, please: ';
	    chomp ($_ = <STDIN>);
	}
	if (/n/i) { die "$outfile not overwritten.\n"; }
	else { open(OUTFILE, ">$outfile"); }
    }
    else { open(OUTFILE, ">$outfile"); }
}

#############
# Calculate #
#############

# Read LDFILE - note that the number of doses should be number of lines - 1

print "What is the number of infected mice per strain per dose (n)? ";
chomp($n = <STDIN>);

@lines = '';
$i = 0;

while (<LDFILE>) {
    chomp($lines[$i++] = $_);
}

# Get doses and survival values

@headers = '';
@headers = split ("\t", $lines[0]);

@survivalmatrix = '';
$j = 0;

while ($j < $#lines) {
    @tempmatrix = split ("\t", $lines[$j+1]);
    $doses[$j] = $tempmatrix[0];
    splice(@tempmatrix, 0, 1);
    push(@survivalmatrix, @tempmatrix);
    $j++;
}

# Transpose matrix

$y = 0;
$z = 0;
@survival = '';
splice (@survivalmatrix, 0, 1);

while ($y < $#headers) {
    $l = 0;
    while ($l <= $#doses) {
	$survival[$z] = $survivalmatrix[$y + ($#headers * $l)];
	$l++;
	$z++;
    }
    $y++;
}

# Tabulate Survival Rates
$k = 0;

while ($k <= $#survival) {
    $m = ($k + 1) % ($#doses + 1);
    if ($m == 0) { $m = $#doses; }
    else { $m--; }

#    print "$k : $m\n";

    $q = 0;
    $totalalive[$k] = 0;
    while ($q <= $m) {
	$totalalive[$k] += $survival[$k - $q];
	$q++;
    }

    $r = 0;
    $totaldead[$k] = 0;
    while ($r <= $#doses - $m) {
	$totaldead[$k] += ($n - $survival[$k + $r]);
	$r++;
    }
    $k++;
}

#foreach $num (@totalalive) {
#    print "$num\n";
#}

#foreach $num (@totaldead) {
#    print "$num\n";
#}

# Calculate Mortality

@mortality = '';
$u = 0;

while ($u <= $#totalalive) {
    $mortality[$u] = $totaldead[$u] / ($totalalive[$u] + $totaldead[$u]);
    $u++;
}

$d = 0;
while ($d <= $#mortality) {
#    print "$totalalive[$d] alive, $totaldead[$d] dead, $mortality[$d] mortality\n";
    $d++;
}

# Calculate LD50

$a = 0;
@ld = '';
$dilution = $doses[0] / $doses[1];
$loval = 0;
$hival = 0;
$propdist = 0;

# sift through mortality values, in groups of # of doses

while ($a < $#headers) {
    $b = 0;
    while ($b < $#doses) {
#	print $b+($a*($#doses+1));
	if ($mortality[$b+($a*($#doses+1))] == 0.5) {
	    $ld[$a] = $doses[$b];
	    last;
	}
	elsif ($mortality[$b+1+($a*($#doses+1))] == 0.5) {
	    $ld[$a] = $doses[$b+1];
	    last;
	}
	elsif ($mortality[$b+($a*($#doses+1))] > 0.5 && $mortality[$b+1+($a*($#doses+1))] < 0.5) {
	    $loval = $mortality[$b+1+($a*($#doses+1))];
	    $hival = $mortality[$b+($a*($#doses+1))];
#	    print "loval $loval hival $hival\n";
	    $propdist = ((0.50-$loval)/($hival-$loval));
#	    print "propdist $propdist\n";
	    $logld = ( (log($doses[$b+1])/log(10) ) + ($propdist * (log($dilution)/log(10)) )  );
	    $ld[$a] =  $dilution ** $logld;
	    last;
	}
	if ($mortality[$b+($a*($#doses+1))] < 0.5) {
	    $ld[$a] = ">$doses[0]";
	}
	if ($mortality[$b+($a*($#doses+1))] > 0.5) {
	    $ld[$a] = "<$doses[$#doses]";
	}
	$b++;
    }
    $a++;
}

$g = 0;
print OUTFILE "Dataset\tLD50\n";
while ($g + 1 <= $#headers) {
    print OUTFILE "$headers[$g + 1]\t$ld[$g]\n";
    $g++;
}

