#!/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++;

}