rsem-control-fdr

   1 #!/usr/bin/perl
   2
   3 use Getopt::Long;
   4 use Pod::Usage;
   5 use strict;
   6
   7 my $hard = 0;
   8 my $soft = 0;
   9 my $help = 0;
  10
  11 GetOptions("hard-threshold" => \$hard,
  12            "soft-threshold" => \$soft,
  13            "h|help" => \$help) or pod2usage(-exitval => 2, -verbose => 2);
  14
  15 pod2usage(-verbose => 2) if ($help == 1);
  16 pod2usage(-msg => "Invalid number of arguments!", -exitval => 2, -verbose => 2) if (scalar(@ARGV) != 3);
  17 pod2usage(-msg => "--hard-threshold and --soft-threshold cannot be set at the same time!", -exitval => 2, -verbose => 2) if ($hard && $soft);
  18
  19 if ($hard == 0 && $soft == 0) { $hard = 1; }
  20
  21 my $fdr = $ARGV[1];
  22
  23 open(INPUT, "$ARGV[0]");
  24 open(OUTPUT, ">$ARGV[2]");
  25
  26 my $header = <INPUT>;
  27 chomp($header);
  28 my @columns = split(/\t/, $header);
  29
  30 my $pos = 0;
  31 while ($pos <= $#columns && $columns[$pos] ne "\"PPDE\"") { ++$pos; }
  32 if ($pos > $#columns) { print "Error: Cannot find column PPDE!\n"; exit(-1); }
  33 ++$pos;
  34
  35 print OUTPUT "$header\n";
  36
  37 my ($n, $sum) = (0, 0);
  38 my $line = "";
  39 while($line = <INPUT>) {
  40     chomp($line);
  41     my @fields = split(/\t/, $line);
  42     my $ppee = 1.0 - $fields[$pos];
  43     if ($hard) {
  44         if ($ppee > $fdr) { last; }
  45         ++$n;
  46         print OUTPUT "$line\n";
  47     }
  48     else {
  49         if ($sum + $ppee > $fdr * ($n + 1)) { last; }
  50         ++$n;
  51         $sum += $ppee;
  52         print OUTPUT "$line\n";
  53     }
  54 }
  55
  56 print "There are $n genes/transcripts reported at FDR = $fdr.\n";
  57
  58 close(INPUT);
  59 close(OUTPUT);
  60
  61 __END__
  62
  63 =head1 NAME
  64
  65 rsem-control-fdr
  66
  67 =head1 SYNOPSIS
  68
  69 rsem-control-fdr [options] input_file fdr_rate output_file
  70
  71 =head1 ARGUMENTS
  72
  73 =over
  74
  75 =item B<input_file>
  76
  77 This should be the result file generated by 'rsem-run-ebseq', which contains all genes/transcripts and their associated statistics.
  78
  79 =item B<fdr_rate>
  80
  81 The desire false discovery rate (FDR).
  82
  83 =item B<output_file>
  84
  85 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.
  86
  87 =back
  88
  89 =head1 OPTIONS
  90
  91 =over
  92
  93 =item B<--hard-threshold>
  94
  95 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)
  96
  97 =item B<--soft-threshold>
  98
  99 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)
 100
 101 =item B<-h/--help>
 102
 103 Show help information.
 104
 105 =back
 106
 107 =head1 DESCRIPTION
 108
 109 This program controls the false discovery rate and reports differentially expressed genes/transcripts.
 110
 111 =head1 EXAMPLES
 112
 113 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':
 114
 115  rsem-control-fdr GeneMat.results 0.05 GeneMat.de.txt
 116
 117 If we instead want to use soft threshold method to obtain more called genes/transcripts, we can use:
 118
 119  rsem-control-fdr --soft-threshold GeneMat.results 0.05 GeneMat.de.txt
 120
 121 =cut