rsem-run-ebseq

   1 #!/usr/bin/env perl
   2
   3 use Getopt::Long;
   4 use Pod::Usage;
   5
   6 use FindBin;
   7 use lib $FindBin::RealBin;
   8 use rsem_perl_utils;
   9
  10 use Env qw(@PATH);
  11 @PATH = ("$FindBin::RealBin/EBSeq", @PATH);
  12
  13 use strict;
  14
  15 my $ngvF = "";
  16 my $help = 0;
  17
  18 GetOptions("ngvector=s" => \$ngvF,
  19            "h|help" => \$help) or pod2usage(-exitval => 2, -verbose => 2);
  20
  21 pod2usage(-verbose => 2) if ($help == 1);
  22 pod2usage(-msg => "Invalid number of arguments!", -exitval => 2, -verbose => 2) if (scalar(@ARGV) != 3);
  23 pod2usage(-msg => "ngvector file cannot be named as #! # is reserved for other purpose!", -exitval => 2, -verbose => 2) if ($ngvF eq "#");
  24
  25 my $command = "";
  26
  27 my @conditions = split(/,/, $ARGV[1]);
  28
  29 pod2usage(-msg => "At least 2 conditions are required for differential expression analysis!", -exitval => 2, -verbose => 2) if (scalar(@conditions) < 2);
  30
  31 if ($ngvF eq "") { $ngvF = "#"; }
  32
  33 $" = " ";
  34 $command = "rsem-for-ebseq-find-DE $FindBin::RealBin/EBSeq $ngvF $ARGV[0] $ARGV[2] @conditions";
  35 &runCommand($command)
  36
  37 __END__
  38
  39 =head1 NAME
  40
  41 rsem-run-ebseq
  42
  43 =head1 SYNOPSIS
  44
  45 rsem-run-ebseq [options] data_matrix_file conditions output_file
  46
  47 =head1 ARGUMENTS
  48
  49 =over
  50
  51 =item B<data_matrix_file>
  52
  53 This file is a m by n matrix. m is the number of genes/transcripts and n is the number of total samples. Each element in the matrix represents the expected count for a particular gene/transcript in a particular sample. Users can use 'rsem-generate-data-matrix' to generate this file from expression result files.
  54
  55 =item B<conditions>
  56
  57 Comma-separated list of values representing the number of replicates for each condition. For example, "3,3" means the data set contains 2 conditions and each condition has 3 replicates. "2,3,3" means the data set contains 3 conditions, with 2, 3, and 3 replicates for each condition respectively.
  58
  59 =item B<output_file>
  60
  61 Output file name.
  62
  63 =back
  64
  65 =head1 OPTIONS
  66
  67 =over
  68
  69 =item B<--ngvector> <file>
  70
  71 This option provides the grouping information required by EBSeq for isoform-level differential expression analysis. The file can be generated by 'rsem-generate-ngvector'. Turning this option on is highly recommended for isoform-level differential expression analysis. (Default: off)
  72
  73 =item B<-h/--help>
  74
  75 Show help information.
  76
  77 =back
  78
  79 =head1 DESCRIPTION
  80
  81 This program is a wrapper over EBSeq. It performs differential expression analysis and can work on two or more conditions. All genes/transcripts and their associated statistcs are reported in one output file. This program does not control false discovery rate and call differential expressed genes/transcripts. Please use 'rsem-control-fdr' to control false discovery rate after this program is finished.
  82
  83 =head1 OUTPUT
  84
  85 =over
  86
  87 =item B<output_file>
  88
  89 This file reports the calculated statistics for all genes/transcripts. It is written as a matrix with row and column names. The row names are the genes'/transcripts' names. The column names are for the reported statistics.
  90
  91 If there are only 2 different conditions among the samples, four statistics (columns) will be reported for each gene/transcript. They are "PPEE", "PPDE", "PostFC" and "RealFC". "PPEE" is the posterior probability (estimated by EBSeq) that a gene/transcript is equally expressed. "PPDE" is the posterior probability that a gene/transcript is differentially expressed. "PostFC" is the posterior fold change (condition 1 over condition2) for a gene/transcript. It is defined as the ratio between posterior mean expression estimates of the gene/transcript for each condition. "RealFC" is the real fold change (condition 1 over condition2) for a gene/transcript.  It is the ratio of the normalized within condition 1 mean count over normalized within condition 2 mean count for the gene/transcript. Fold changes are calculated using EBSeq's 'PostFC' function. The genes/transcripts are reported in descending order of their "PPDE" values.
  92
  93 If there are more than 2 different conditions among the samples, the output format is different. For differential expression analysis with more than 2 conditions, EBSeq will enumerate all possible expression patterns (on which conditions are equally expressed and which conditions are not). Suppose there are k different patterns, the first k columns of the output file give the posterior probability of each expression pattern is true. Patterns are defined in a separate file, 'output_file.pattern'. The k+1 column gives the maximum a posteriori (MAP) expression pattern for each gene/transcript. The k+2 column gives the posterior probability that not all conditions are equally expressed (column name "PPDE"). The genes/transcripts are reported in descending order of their "PPDE" column values. For details on how EBSeq works for more than 2 conditions, please refer to EBSeq's manual.
  94
  95 =item B<output_file.pattern>
  96
  97 This file is only generated when there are more than 2 conditions. It defines all possible expression patterns over the conditions using a matrix with names. Each row of the matrix refers to a different expression pattern and each column gives the expression status of a different condition. Two conditions are equally expressed if and only if their statuses are the same.
  98
  99 =item B<output_file.condmeans>
 100
 101 This file is only generated when there are more than 2 conditions. It gives the normalized mean count value for each gene/transcript at each condition. It is formatted as a matrix with names. Each row represents a gene/transcript and each column represent a condition. The order of genes/transcripts is the same as 'output_file'. This file can be used to calculate fold changes between conditions which users are interested in.
 102
 103 =back
 104
 105 =head1 EXAMPLES
 106
 107 1) We're interested in isoform-level differential expression analysis and there are two conditions. Each condition has 5 replicates. We have already collected the data matrix as 'IsoMat' and generated ngvector as 'ngvector.ngvec':
 108
 109  rsem-run-ebseq --ngvector ngvector.ngvec IsoMat 5,5 IsoMat.results
 110
 111 The results will be in 'IsoMat.results'.
 112
 113 2) We're interested in gene-level analysis and there are 3 conditions. The first condition has 3 replicates and the other two has 4 replicates each. The data matrix is named as 'GeneMat':
 114
 115  rsem-run-ebseq GeneMat 3,4,4 GeneMat.results
 116
 117 Three files, 'GeneMat.results', 'GeneMat.results.pattern', and 'GeneMat.results.condmeans', will be generated.
 118
 119 =cut
 120
 121
 122
 123
 124