#!/usr/bin/perl -w # DO NOT EDIT the /cluster/bin/scripts copy of this file -- # edit ~/kent/src/utils/lodToBedScore instead. # $Id: lodToBedScore,v 1.4 2004/10/11 00:18:03 acs Exp $ use strict; # This script does Adam's log transform of input scores into the bed score # "spectrum" for browser display: 0-1000. The transform is # f(x) = a * log x + b, s.t. f(x_med) = 300 and f(x_max) = 1000. # where x_max is the max input score and x_med is the median input score. # Solving for a and b, you obtain # b = (300 log x_max - 1000 log x_med) / (log x_max - log x_med) # a = (1000 - b) / log x_max. sub usage { die " Usage: $0 [input.bed] Transform scores in input.bed (which must have at least 5 columns, whitespace-delimited) by the function f(x) = a * log x + b, where f(x_med) = 300 and f(x_max) = 1000 to stdout. Inputs must be positive, and max must be greater than median, in order for this transform to work. \n"; } if (grep /^-h/, @ARGV) { &usage(); } # slurp input into memory (duplicate scores so they can be sorted separately): my @lines = (); my @scores = (); while (<>) { next if (/^#/); chomp; my @words = split(); if (scalar(@words) < 5) { die "\nError: input must have at least 5 columns. (input line $.)\n\n"; } if ($words[4] <= 0) { printf STDERR "\nWarning: original scores must be positive; value of %f will be transformed to 0 (input line $.)\n\n", $words[4]; } else { # don't push non-positive scores push @scores, $words[4]; } push @lines, \@words; } # sort scores; extract and check max and median. @scores = sort {$a <=> $b} @scores; my $n = scalar(@scores); my $max = $scores[$n-1]; my $med = $scores[int($n / 2)]; if ($max <= $med) { die " Error: max score ($max) must be greater than median ($med) in order for this transform to apply. \n"; } # calculate coefficients a and b. my $b = ( ((300 * log($max)) - (1000 * log($med))) / (log($max) - log($med)) ); my $a = (1000 - $b) / log($max); # Do the transform on the score column of saved input: foreach my $wRef (@lines) { my @words = @{$wRef}; if ($words[4] > 0) { $words[4] = int(($a * log($words[4])) + $b); } if ($words[4] < 0) { $words[4] = 0 }; # applies either if $words[4] was # originally negative or if it is made # negative by the transformation print join("\t", @words) . "\n"; }