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

# Example 3: ranklods.pl
#	Rank individual family LOD scores
#
# 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 15, 1999
#
# Usage: ranklods.pl 
#	[marker number]
#	[LINKAGE/FASTLINK-MLINK output, usually called final.out]
#	[post-MAKEPED pedigree file]
#

$usage ="Usage: ranklods.pl [marker] [final.out] [.ped]";

if($#ARGV != 2 ){ print "$usage\n"; exit; }

# marker locus number in the order found in in parameter & pedigree file (
# assuming disease locus is the first 
$locus = $ARGV[0];  

# file name for the linkage output  files (e.g. final.out)
$outFile = $ARGV[1];

# file name for post-makeped pedigree file
$pedFile = $ARGV[2];

open(OUT, "< $outFile") || die "Can't open $outFile file to read:$!\n";
open(PED, "< $pedFile") || die "Can't open $pedFile file to read:$!\n";


# find "Locus Order" line for marker locus of interst
while(<OUT>) {
    chomp;
    last if(m/(Locus Order) [^:]+: (1 $locus)/);
}


while(<OUT>) {
    # exit loop at next "Locus Order" line
    last if(m/(Locus Order) [^:]+:/); 
    s/^\s+//g;
    if(m/THETAS/) {
	$theta = (split)[1];
	next;
    }
    next unless(m/LOD=/); # feature of FASTLINK output only
    # $key => family ID and $lod => family lod score
    ($key, $lod) = (split)[0,4];

    # for each family, store max lod sofar in a table
    if(!defined $lodHash{$key}) {
	$lodHash{$key} = {
	    ID => $key,
	    THETA => $theta,
	    LOD => $lod
	    }
    }
    elsif($lodHash{$key}->{LOD} < $lod) {
	$lodHash{$key}->{THETA} = $theta;
	$lodHash{$key}->{LOD} = $lod;
    }
}

# store pedigree file into array
my @ped = grep (/\n/, <PED>);
$prevPostID =  -1;

# match pedigree ID from lod table (lodHash) with pre-makeped pedigree ID 
# assumed to be at end of each line in the pedigree file ($pedFile)

foreach $i (0..$#ped) {
    $ped[$i] =~ s/^\s+//; # remove space in begining of line
    @pedLine = split(/Ped: /, $ped[$i]);
    $preID = (split(/\s+/, $pedLine[1]))[0]; # pre-makeped pedigree ID
    $postID = (split(/\s+/, $pedLine[0]))[0]; # post-makeped pedigree ID
    if($postID ne $prevPostID) {
	$preToPost{$postID} = $preID;
	$prevPostID = $postID;
    }
}

# sort  and print the per family maximum lod score in decending order
@sortByLodRef = sort {$$b{LOD} <=> $$a{LOD} } values %lodHash;
print "Locus $locus per family ranked maximum lod scores:\n";
for $fam (@sortByLodRef) {
    $preID = $$fam{ID};
    print "$preID <=> $preToPost{$preID}\t\t",
    "lod = $$fam{LOD}\t\t",
    "theta = $$fam{THETA}\n";
}