#!/usr/bin/perl
# AUTHOR: Timothy L. Bailey
# CREATE DATE: May 7, 2010
#use strict;
use File::Basename;

$PGM = $0;			# name of program
$PGM =~ s#.*/##;                # remove part up to last slash
#@args = @ARGV;			# arguments to program
$| = 1;				# flush after all prints
$SIG{'INT'} = \&cleanup;	# interrupt handler
# Note: so that interrupts work, always use for system calls:
# 	if ($status = system("$command")) {&cleanup($status)}

# requires
push(@INC, split(":", $ENV{'PATH'}));	# look in entire path

# defaults
my $seed = 1;
my $bootstraps = 1000;

my $usage = <<USAGE;		# usage message
  USAGE:
	$PGM [options]

	[-seed <seed>]		seed for random numbers; default: $seed
	[-bootstraps <boot>]	number of bootstraps to perform 
				while computing pi_0; default: 1000

	Read an AMA output file and output each sequence name, p-value
	and q-value.  At the end of file, the value of pi_0 (number of
	sequences not showing significant binding to motif) is shown,
	unless there are fewer than 100 p-values in the input.

	Reads standard input.
	Writes standard output.

        Copyright
        (2010) The University of Queensland
        All Rights Reserved.
        Author: Timothy L. Bailey
USAGE

$nargs = 0;			# number of required args
if ($#ARGV+1 < $nargs) { &print_usage("$usage", 1); }

# get input arguments
while ($#ARGV >= 0) {
  $_ = shift;
  if ($_ eq "-h") {				# help
    &print_usage("$usage", 0);
  } elsif ($_ eq "-seed") {
    $seed = shift;
  } elsif ($_ eq "-bootstraps") {
    $bootstraps = shift;
  } else {
    &print_usage("$usage", 1);
  }
}

# open a pipe to the qvalue program, reading standard input
$pzfile = "$PGM.$$.pzfile.tmp";
open(QVALUE, "|qvalue --append --column 2 --pi-zero-file $pzfile --seed $seed --bootstraps $bootstraps -") || 
  die("Can't open qvalue program.\n");

# compute the qvalues
print "# sequence_name\t\t\tp-value\t\tq-value\n";
my($pvalue, @w, $name);
my $npvalues = 0;
while (<STDIN>) {
  if (/pvalue="([^"]*)"/) {
    $pvalue = $1; 
    /name="([^"]*)"/; 
    $name = $1; 
    print(QVALUE "$name\t$pvalue\n");
    $npvalues++;
  }
};
close(QVALUE);

printf("# Fraction of sequences not showing significant binding: ");
system("cat $pzfile") if ($npvalues >= 100);
unlink $pzfile;

$status = 1;

# cleanup files
&cleanup($status, "");
 
################################################################################
#                       Subroutines                                            #
################################################################################
 
################################################################################
#
#       print_usage
#
#	Print the usage message and exit.
#
################################################################################
sub print_usage {
  my($usage, $status) = @_;
 
  if (-c STDOUT) {			# standard output is a terminal
    open(C, "| more");
    print C $usage;
    close C;
  } else {				# standard output not a terminal
    print STDERR $usage;
  }

  exit $status;
}
 
################################################################################
#       cleanup
#
#       cleanup stuff
#
################################################################################
sub cleanup {
  my($status, $msg) = @_;
  if ($status eq "INT") {
    exit(1);
  } else {
    if ($status && "$msg") {print STDERR "$msg: $status\n";}
  }
}