Updated boost to v1.55.0

[rsem.git] / rsem-calculate-expression

diff --git a/rsem-calculate-expression b/rsem-calculate-expression
index 03811afb0c8b79c1ae348ce01883820908e0c4b8..ccfc643e1c0cbb2963fcc34159193dd72cd27103 100755 (executable)
--- a/rsem-calculate-expression
+++ b/rsem-calculate-expression
@@ -57,7 +57,6 @@ my $genGenomeBamF = 0;
 my $sampling = 0;
 my $calcPME = 0;
 my $calcCI = 0;
-my $var_opt = 0; # temporarily, only for internal use
 my $quiet = 0;
 my $help = 0;
 
@@ -73,6 +72,8 @@ my $bowtie2_mismatch_rate = 0.1;
 my $bowtie2_k = 200;
 my $bowtie2_sensitivity_level = "sensitive"; # must be one of "very_fast", "fast", "sensitive", "very_sensitive"
 
+my $seed = "NULL";
+
 my $version = 0;
 
 my $mTime = 0;
@@ -122,10 +123,10 @@ GetOptions("keep-intermediate-files" => \$keep_intermediate_files,
           "output-genome-bam" => \$genGenomeBamF,
           "sampling-for-bam" => \$sampling,
           "calc-pme" => \$calcPME,
-          "var" => \$var_opt,
           "calc-ci" => \$calcCI,
           "ci-memory=i" => \$NMB,
           "samtools-sort-mem=s" => \$SortMem,
+          "seed=i" => \$seed,
           "time" => \$mTime,
           "version" => \$version,
           "q|quiet" => \$quiet,
@@ -160,6 +161,7 @@ pod2usage(-msg => "Number of threads should be at least 1!\n", -exitval => 2, -v
 pod2usage(-msg => "Seed length should be at least 5!\n", -exitval => 2, -verbose => 2) if ($L < 5);
 pod2usage(-msg => "--sampling-for-bam cannot be specified if --no-bam-output is specified!\n", -exitval => 2, -verbose => 2) if ($sampling && !$genBamF);
 pod2usage(-msg => "--output-genome-bam cannot be specified if --no-bam-output is specified!\n", -exitval => 2, -verbose => 2) if ($genGenomeBamF && !$genBamF);
+pod2usage(-msg => "The seed for random number generator must be a non-negative 32bit integer!\n", -exitval => 2, -verbose => 2) if (($seed ne "NULL") && ($seed < 0 || $seed > 0xffffffff));
 
 if ($L < 25) { print "Warning: the seed length set is less than 25! This is only allowed if the references are not added poly(A) tails.\n"; }
 
@@ -335,14 +337,23 @@ print OUTPUT "$mean $sd\n";
 print OUTPUT "$L\n";
 close(OUTPUT);  
 
+my @seeds = ();
+if ($seed ne "NULL") { 
+    srand($seed); 
+    for (my $i = 0; $i < 3; $i++) {
+       push(@seeds, int(rand(1 << 32))); 
+    }
+}
+
 $command = "rsem-run-em $refName $read_type $sampleName $imdName $statName -p $nThreads";
 if ($genBamF) { 
     $command .= " -b $samInpType $inpF";
     if ($fn_list ne "") { $command .= " 1 $fn_list"; }
     else { $command .= " 0"; }
     if ($sampling) { $command .= " --sampling"; }
+    if ($seed ne "NULL") { $command .= " --seed $seeds[0]"; }
 }
-if ($calcPME || $var_opt || $calcCI) { $command .= " --gibbs-out"; }
+if ($calcPME || $calcCI) { $command .= " --gibbs-out"; }
 if ($quiet) { $command .= " -q"; }
 
 &runCommand($command);
@@ -377,10 +388,10 @@ if ($mTime) { $time_end = time(); $time_rsem = $time_end - $time_start; }
 
 if ($mTime) { $time_start = time(); }
 
-if ($calcPME || $var_opt || $calcCI ) {
+if ($calcPME || $calcCI ) {
     $command = "rsem-run-gibbs $refName $imdName $statName $BURNIN $NCV $SAMPLEGAP";
     $command .= " -p $nThreads";
-    if ($var_opt) { $command .= " --var"; }
+    if ($seed ne "NULL") { $command .= " --seed $seeds[1]"; }
     if ($quiet) { $command .= " -q"; }
     &runCommand($command);
 }
@@ -405,6 +416,7 @@ if ($calcPME || $calcCI) {
 if ($calcCI) {
     $command = "rsem-calculate-credibility-intervals $refName $imdName $statName $CONFIDENCE $NCV $NSPC $NMB";
     $command .= " -p $nThreads";
+    if ($seed ne "NULL") { $command .= " --seed $seeds[2]"; }
     if ($quiet) { $command .= " -q"; }
     &runCommand($command);
 
@@ -526,6 +538,10 @@ Generate a BAM file, 'sample_name.genome.bam', with alignments mapped to genomic
 
 When RSEM generates a BAM file, instead of outputing all alignments a read has with their posterior probabilities, one alignment is sampled according to the posterior probabilities. The sampling procedure includes the alignment to the "noise" transcript, which does not appear in the BAM file. Only the sampled alignment has a weight of 1. All other alignments have weight 0. If the "noise" transcript is sampled, all alignments appeared in the BAM file should have weight 0. (Default: off)
 
+=item B<--seed> <uint32>
+
+Set the seed for the random number generators used in calculating posterior mean estimates and credibility intervals. The seed must be a non-negative 32 bit interger. (Default: off)
+
 =item B<--calc-pme>
 
 Run RSEM's collapsed Gibbs sampler to calculate posterior mean estimates. (Default: off) 
@@ -688,7 +704,7 @@ File containing isoform level expression estimates. The first line
 contains column names separated by the tab character. The format of
 each line in the rest of this file is:
 
-transcript_id gene_id length effective_length expected_count TPM FPKM IsoPct [pme_expected_count pme_TPM pme_FPKM IsoPct_from_pme_TPM TPM_ci_lower_bound TPM_ci_upper_bound FPKM_ci_lower_bound FPKM_ci_upper_bound]
+transcript_id gene_id length effective_length expected_count TPM FPKM IsoPct [posterior_mean_count posterior_standard_deviation_of_count pme_TPM pme_FPKM IsoPct_from_pme_TPM TPM_ci_lower_bound TPM_ci_upper_bound FPKM_ci_lower_bound FPKM_ci_upper_bound]
 
 Fields are separated by the tab character. Fields within "[]" are
 optional. They will not be presented if neither '--calc-pme' nor
@@ -734,9 +750,11 @@ transcript's abandunce over its parent gene's abandunce. If its parent
 gene has only one isoform or the gene information is not provided,
 this field will be set to 100.
 
-'pme_expected_count', 'pme_TPM', 'pme_FPKM' are posterior mean
-estimates calculated by RSEM's Gibbs sampler. 'IsoPct_from_pme_TPM' is
-the isoform percentage calculated from 'pme_TPM' values.
+'posterior_mean_count', 'pme_TPM', 'pme_FPKM' are posterior mean
+estimates calculated by RSEM's Gibbs
+sampler. 'posterior_standard_deviation_of_count' is the posterior
+standard deviation of counts. 'IsoPct_from_pme_TPM' is the isoform
+percentage calculated from 'pme_TPM' values.
 
 'TPM_ci_lower_bound', 'TPM_ci_upper_bound', 'FPKM_ci_lower_bound' and
 'FPKM_ci_upper_bound' are lower(l) and upper(u) bounds of 95%
@@ -749,7 +767,7 @@ File containing gene level expression estimates. The first line
 contains column names separated by the tab character. The format of
 each line in the rest of this file is:
 
-gene_id transcript_id(s) length effective_length expected_count TPM FPKM [pme_expected_count pme_TPM pme_FPKM TPM_ci_lower_bound TPM_ci_upper_bound FPKM_ci_lower_bound FPKM_ci_upper_bound]
+gene_id transcript_id(s) length effective_length expected_count TPM FPKM [posterior_mean_count posterior_standard_deviation_of_count pme_TPM pme_FPKM TPM_ci_lower_bound TPM_ci_upper_bound FPKM_ci_lower_bound FPKM_ci_upper_bound]
 
 Fields are separated by the tab character. Fields within "[]" are
 optional. They will not be presented if neither '--calc-pme' nor
@@ -774,7 +792,7 @@ allele-specific expression calculation. The first line
 contains column names separated by the tab character. The format of
 each line in the rest of this file is:
 
-allele_id transcript_id gene_id length effective_length expected_count TPM FPKM AlleleIsoPct AlleleGenePct [pme_expected_count pme_TPM pme_FPKM AlleleIsoPct_from_pme_TPM AlleleGenePct_from_pme_TPM TPM_ci_lower_bound TPM_ci_upper_bound FPKM_ci_lower_bound FPKM_ci_upper_bound]
+allele_id transcript_id gene_id length effective_length expected_count TPM FPKM AlleleIsoPct AlleleGenePct [posterior_mean_count posterior_standard_deviation_of_count pme_TPM pme_FPKM AlleleIsoPct_from_pme_TPM AlleleGenePct_from_pme_TPM TPM_ci_lower_bound TPM_ci_upper_bound FPKM_ci_lower_bound FPKM_ci_upper_bound]
 
 Fields are separated by the tab character. Fields within "[]" are
 optional. They will not be presented if neither '--calc-pme' nor