Added support for DE analysis on multiple conditions via running EBSeq

[rsem.git] / rsem-control-fdr

diff --git a/rsem-control-fdr b/rsem-control-fdr
new file mode 100755 (executable)
index 0000000..86cfb93
--- /dev/null
+++ b/rsem-control-fdr
@@ -0,0 +1,121 @@
+#!/usr/bin/perl
+
+use Getopt::Long;
+use Pod::Usage;
+use strict;
+
+my $hard = 0;
+my $soft = 0;
+my $help = 0;
+
+GetOptions("hard-threshold" => \$hard,
+          "soft-threshold" => \$soft,
+          "h|help" => \$help) or pod2usage(-exitval => 2, -verbose => 2);
+
+pod2usage(-verbose => 2) if ($help == 1);
+pod2usage(-msg => "Invalid number of arguments!", -exitval => 2, -verbose => 2) if (scalar(@ARGV) != 3);
+pod2usage(-msg => "--hard-threshold and --soft-threshold cannot be set at the same time!", -exitval => 2, -verbose => 2) if ($hard && $soft);
+
+if ($hard == 0 && $soft == 0) { $hard = 1; }
+
+my $fdr = $ARGV[1];
+
+open(INPUT, "$ARGV[0]");
+open(OUTPUT, ">$ARGV[2]");
+
+my $header = <INPUT>;
+chomp($header);
+my @columns = split(/\t/, $header);
+
+my $pos = 0;
+while ($pos <= $#columns && $columns[$pos] ne "\"PPDE\"") { ++$pos; }
+if ($pos > $#columns) { print "Error: Cannot find column PPDE!\n"; exit(-1); }
+++$pos;
+
+print OUTPUT "$header\n";
+
+my ($n, $sum) = (0, 0);
+my $line = "";
+while($line = <INPUT>) {
+    chomp($line);
+    my @fields = split(/\t/, $line);
+    my $ppee = 1.0 - $fields[$pos];
+    if ($hard) {
+       if ($ppee > $fdr) { last; }
+       ++$n;
+       print OUTPUT "$line\n";
+    }
+    else {
+       if ($sum + $ppee > $fdr * ($n + 1)) { last; }
+       ++$n;
+       $sum += $ppee;
+       print OUTPUT "$line\n";
+    }
+}
+
+print "There are $n genes/transcripts reported at FDR = $fdr.\n";
+
+close(INPUT);
+close(OUTPUT);
+
+__END__
+
+=head1 NAME
+
+rsem-control-fdr
+
+=head1 SYNOPSIS
+
+rsem-control-fdr [options] input_file fdr_rate output_file
+
+=head1 ARGUMENTS
+
+=over
+
+=item B<input_file>
+
+This should be the result file generated by 'rsem-run-ebseq', which contains all genes/transcripts and their associated statistics.
+
+=item B<fdr_rate>
+
+The desire false discovery rate (FDR). 
+
+=item B<output_file>
+
+This file is a subset of the 'input_file'. It only contains the genes/transcripts called as differentially expressed (DE). When more than 2 conditions exist, DE is defined as not all conditions are equally expressed.
+
+=back
+
+=head1 OPTIONS
+
+=over
+
+=item B<--hard-threshold> 
+
+Use hard threshold method to control FDR. If this option is set, only those genes/transcripts with their PPDE >= 1 - fdr_rate are called as DE. (Default: on)
+
+=item B<--soft-threshold>
+
+Use soft threshold method to control FDR. If this option is set, this program will try to report as many genes/transcripts as possible, as long as their average PPDE >= 1 - fdr_rate. (Default: off)
+
+=item B<-h/--help> 
+
+Show help information.
+
+=back
+
+=head1 DESCRIPTION
+
+This program controls the false discovery rate and reports differentially expressed genes/transcripts. 
+
+=head1 EXAMPLES
+
+We assume that we have 'GeneMat.results' as input. We want to control FDR at 0.05 using hard threshold method and name the output file as 'GeneMat.de.txt':
+
+ rsem-control-fdr GeneMat.results 0.05 GeneMat.de.txt
+
+If we instead want to use soft threshold method to obtain more called genes/transcripts, we can use:
+
+ rsem-control-fdr --soft-threshold GeneMat.results 0.05 GeneMat.de.txt
+
+=cut