#!/usr/local/bin/perl -w

# Example 2: lia_panic.pl
#	Creating a liability class column in LINKAGE pedigree
#	file as well as assigning the corresponding penetrance 
#	tables in the LINKAGE parameter file.
#	Case 2: Panic disorder.
#
# Ref: Li and Haghighi, "Perl as a tool for linkage analysis",
# Am J Hum Genet, suppl,65: A260 (1999)
#
# Contact: perl@linkage.rockefeller.edu
#
# Last update: October 17, 1999
#
# Two Usages: 
#	lia_panic.pl -d  
#        [class definition file] 
#	 [LINKAGE pre-MAKEPED pedigree file]
#
#	lia_panic.pl -p
#	 [class definition file]
#	 [LINKAGE pre-MAKEPED pedigree file]
#	 [dianostic code file]
#	 [-nl (no liability class column) or -wl (with liability column)]

$usage ="Usage: lia.pl -d [class file] [.ped]\nUsage: lia.pl -p [class file] [.ped file] [diag file] [-nl/-wl]";


# parse command line args
$option = $ARGV[0];
$paraFile = $ARGV[1]; 
&usage("Unknown option <$ARGV[0] >") unless($option =~ m/^-/);


$female = 2;
$male = 1;
$unisex = -1;
# $noAgeInfo = 0.0;

# parse command line args
if($ARGV[0] =~ m/^-/) {
    if($ARGV[0] eq "-d") { # recode LINKAGE datafile only
	$paraFile = $ARGV[1];
	$dataFile = $ARGV[2];
	open(IN, "< $paraFile") 
	    || die "Can't open file $paraFile to read:$!\n"; 
	open(DF, "< $dataFile") 
	    || die "Can't open file $dataFile to read:$!\n"; 

	@paraArray = grep(/\n/, <IN>);
	@dataArray = grep(/\n/, <DF>);

	chop($paraArray[3]);
	$numLiab = (split(/=/, $paraArray[3]))[1];
	%liabHash = &initLiab();
	

	print "@dataArray[0..3]";
	$dataArray[5] =~ s/^\s+//;
	$oldNumLiab = (split(/\s+/, $dataArray[5]))[0];
	$startMrk = $oldNumLiab + 6;
	

	$alleleIndx = $numLiab + 4;
	$paraArray[$alleleIndx] =~ s/^\s+//;
	$dFrq = (split(/\s+/, $paraArray[$alleleIndx]))[0];
	$nFrq = 1 - $dFrq;
	print "$nFrq $dFrq\n";
	print "$numLiab \n";
	foreach $key (1..$numLiab) {
	    print "$liab{$key}->{PENE}\n";
	}
	print "@dataArray[$startMrk..$#dataArray]";

    }
    elsif($ARGV[0] eq "-p") { # recode LINKAGE pedfile
	$paraFile = $ARGV[1];
	$pedFile = $ARGV[2];
	$codeFile = $ARGV[3];
	open(CF, "< $codeFile") 
	    || die "Can't open file $codeFile to read:$!\n";
	open(PF, "< $pedFile") 
	    || die "Can't open file $pedFile to read:$!\n";
	open(IN, "< $paraFile") 
	    || die "Can't open file $paraFile to read:$!\n";
	@codeArray = grep(/\n/, <CF>);
	# print @codeArray;
	@pedArray = grep(/\n/, <PF>);
	# print @pedArray;
	@paraArray = grep(/\n/, <IN>);
	# print @paraArray;
	chop($paraArray[3]);
	$numLiab = (split(/=/, $paraArray[3]))[1];
	%liabHash = &initLiab();
	%codeHash = &newLiab();
	foreach $p (@pedArray) {
	    @person = split(/\s+/, $p);
	    $pkey = $person[0]."|".$person[1];
	    if(defined $codeHash{$pkey}) {
		if($numLiab == 1) {
		    if($ARGV[4] eq "-wl") {
			print "@person[0..4]",
			" $codeHash{$pkey}->{AFFSTATUS}",
			" @person[7..$#person]\n";
		    }
		    elsif($ARGV[4] eq "-nl") {
			print "@person[0..4]",
			" $codeHash{$pkey}->{AFFSTATUS}",
			" @person[6..$#person]\n";
		    }
		    else {
			&usage("Incorrect pedigree LIAB option"); 
		    }
		}
		else {
                    if($ARGV[4] eq "-wl") {
			print "@person[0..4] $codeHash{$pkey}->{AFFSTATUS}",
			" $codeHash{$pkey}->{LIAB} @person[7..$#person]\n";
                    }
                    elsif($ARGV[4] eq "-nl") {
			print "@person[0..4] $codeHash{$pkey}->{AFFSTATUS}",
			" $codeHash{$pkey}->{LIAB} @person[6..$#person]\n";
                    }
                    else {
                        &usage("Incorrect pedigree LIAB option");
                    }
		}
		next;
	    }
	    else {
		print STDERR "ERROR:  There is no CODE for",
		" person $person[1] in pedigree $person[0]. \n";
		# exit(-1);
	    } 
	} # end foreach $p (@pedArray)
    } # end elsif "-p"
    else { # unknown option
	&usage("The option ($ARGV[0]) is unknown!");
    }
}
else {
    &usage("Incorrect command line option (use -d or -p)");
}




sub usage {
    my $errString = shift;
    print STDERR "$errString\n";
    print STDERR "This program uses the following command line formats:\n";
    print STDERR "\t1) Recode LINKAGE datafile only => -d \n",
    "\t<paramater file> <datafile>\n";
    print STDERR "\t2) Recode LINKAGE pedfile only => -p \n",
    "\t<paramater file> <pedfile> <code file>\n";
    print STDERR "\t3) Recode LINKAGE datafile & pedfile => \n",
    "\t<paramater file> <datafile> <pedfile> <codefile>\n";
    print STDERR "NOTE-1: All pedigree files must be in\n",
    "\t pre- MAKEPED LINKAGE format.\n";
    print STDERR "NOTE-2: All penetrances in the paramater\n",
    "\t file  must be in the order ++, D+ , DD.\n";
    exit(-1);
}

sub initLiab {
    foreach $i (1..$numLiab) {
	$paraIndex = $i+3;
	chop($paraArray[$paraIndex]);
	local @l = split(/,/, $paraArray[$paraIndex]);
	$sex = $female if($l[1] =~ m/F|f/);
	$sex = $male if($l[1] =~ m/M|m/);
	$sex = $unisex if($l[1] =~ m/U|u/);
	$l[2] =~ s/^\s+// ;
	@ageRange = split(/\s+/, $l[2]);
	$liab{$i} = {
	    AFFSTATUS => $l[0],
	    SEX       => $sex,
	    AGELOW  => $ageRange[0],
	    AGEUP  =>  $ageRange[1],
	    PENE      => $l[3]
	    }
    }
    return(%liab);
}


sub newLiab {
    foreach $i (0..2) {
	$stat = (split(/=/, $paraArray[$i]))[1];
	$stat =~ s/^\s+//;
	@aff = split(/\s+/, $stat) 
	    if($paraArray[$i] =~ m/^AFFECTED|affected/);
	@unaff = split(/\s+/, $stat) 
	    if($paraArray[$i] =~ m/^UNAFFECTED|unaffected/);
	@unk = split(/\s+/, $stat) 
	    if($paraArray[$i] =~ m/^UNKNOWN|unknown/);
    }
    foreach $pCode (@codeArray) {
        last if($pCode =~ /^\s+$/);
	@person = split(/\s+/,$pCode);
	$affStatus = &findAffStatus(\@aff, 
				    \@unaff, 
				    \@unk, 
				    $person[4]);
	if($affStatus == -1) {
	    print STDERR "ERROR:  The diagnostic code for the following",
	    " person is unknown: $pCode\n";
	}
	
	$refaff_liab = &getCodeH($affStatus, \@unaff, \@unk, \@person);
	$key = $person[0]."|".$person[1];

	if(defined $codeH{$key}) {
	    print STDERR "ERROR: Person $person[1] in pedigree $person[0] has",
	    "multiple records in the code file.\n";
	    exit(-1);
	}
	else {
	    $codeH{$key} = {
		AFFSTATUS => @$refaff_liab[0],
		LIAB => @$refaff_liab[1]
		}
	}
    }
    return(%codeH);
}

sub findAffStatus {
    my ($raffList, $runaffList, $runkList, $diag) = @_;
	return(2) if(grep /^$diag$/, @$raffList);
	return(1) if(grep /^$diag$/, @$runaffList);
	return(0) if(grep /^$diag$/, @$runkList);
	return(-1); # when not diagnoses match
}


sub findLiab {
    my ($affFlag, $refp) = @_;
    local $diag = @$refp[4];
    my ($key, @keys);
    @keys = sort {$a <=> $b} keys %liabHash;
    foreach $key (@keys) {
	$rliab = $liabHash{$key};
	if($rliab->{AFFSTATUS}=~ m/$diag/)  {
	    if((@$refp[2] == $female) || (@$refp[2] == $male)) {
		if(($rliab->{SEX} == @$refp[2]) || 
		   ($rliab->{SEX} == $unisex)){
		    # call age func.
		    my $ageStatus = &checkAge($affFlag, 
					      $rliab->{AGELOW}, 
					      $rliab->{AGEUP}, 
					      $refp);
		    return($key)  if($ageStatus == 1);
		    next if($ageStatus == 0);
		    last; # not liab. class matched 
		}
	    }
	    else { # person is not male or female
		print STDERR "ERROR:  The following person has incorrect",
		" sex:\nPERSON => @$refp\n";
		exit(-1);
	    }
	}
    }
    return($numLiab + 1);
}

sub checkAge {
    my ($affFlag, $low, $up, $refp) = @_;
    if($affFlag == 2) {
	if(@$refp[5] > 0) {
	     return(1) if((@$refp[5] >= $low) && (@$refp[5] <= $up)); 
	     return(0); 
	}
	elsif(@$refp[6] > 0) {
	    return(1) if((@$refp[6] >= $low) && (@$refp[6] <= $up));
	    return(0); 
	}
	else { # no age of onset known 
            #cover cases where age not important (e.g. 0 to ?)
	    return(1) if($low <= 0); 
	    return(-1); 
	} 
    }
    else { # person not affected
	$interviewAge = sprintf("%.0f", @$refp[3]);
	return(1) if(($interviewAge >= $low) && ($interviewAge <= $up));
	return(0); 
    }
}


sub getCodeH {
    my ($affStatus, $refunaff, $refunk, $refp) = @_;
    $liab = &findLiab($affStatus, $refp);
    while(1) {
	if($liab > $numLiab) {
	    if($numLiab == 1) {
		if($affStatus == 2) { 
		    $affStatus = 1;
		    print STDERR "WARNING:  Coded as Unaffected: @$refp\n";
		}
		last;
	    }
	    if($affStatus == 2) { # if person is affected
		print STDERR "WARNING:  Coded as Unaffected: @$refp\n";
		$affStatus = 1;
		@$refp[4] = @$refunaff[0];
		&getCodeH($affStatus, $refunaff, $refunk, $refp);
	    }
	    elsif($affStatus == 1) { # if person is unaffected
		print STDERR "WARNING:  Coded as Unknown: @$refp\n";
		$affStatus = 0;
		@$refp[4] = @$refunk[0];
		&getCodeH($affStatus, $refunaff, $refunk, $refp);
	    }
	    else { # if person has diag. unknown
		last;
	    }
	} # end if $liab > $numLiab
	else {
	    last;
	}
    } # end while
    $aff_liab[0] = $affStatus;
    $aff_liab[1] = $liab;
    return(\@aff_liab);
}