add optional field XS for cufflinks

[rsem.git] / rsem-calculate-expression

diff --git a/rsem-calculate-expression b/rsem-calculate-expression
index 93fabce4686c8d0b7e4fafeae279c6262e7d412e..14679fff7340ab459ccbe8def6a61139cf9e548c 100755 (executable)
--- a/rsem-calculate-expression
+++ b/rsem-calculate-expression
@@ -45,6 +45,7 @@ my $B = 20;
 
 my $nThreads = 1;
 my $genBamF = 0;
+my $sampling = 0;
 my $calcCI = 0;
 my $quiet = 0;
 my $help = 0;
@@ -83,6 +84,7 @@ GetOptions("keep-intermediate-files" => \$keep_intermediate_files,
           "num-rspd-bins=i" => \$B,
           "p|num-threads=i" => \$nThreads,
           "out-bam" => \$genBamF,
+          "sampling-for-bam" => \$sampling,
           "calc-ci" => \$calcCI,
           "ci-memory=i" => \$NMB,
           "time" => \$mTime,
@@ -110,6 +112,8 @@ pod2usage(-msg => "Min fragment length should be at least 1!", -exitval => 2, -v
 pod2usage(-msg => "Min fragment length should be smaller or equal to max fragment length!", -exitval => 2, -verbose => 2) if ($minL > $maxL);
 pod2usage(-msg => "The memory allocated for calculating credibility intervals should be at least 1 MB!\n", -exitval => 2, -verbose => 2) if ($NMB < 1);
 pod2usage(-msg => "Number of threads should be at least 1!\n", -exitval => 2, -verbose => 2) if ($nThreads < 1);
+pod2usage(-msg => "Seed length should be at least 25!\n", -exitval => 2, -verbose => 2) if ($L < 25);
+pod2usage(-msg => "--sampling-for-bam cannot be specified if --out-bam is not specified!\n", -exitval => 2, -verbose => 2) if ($sampling && !$genBamF);
 
 if ($strand_specific) { $probF = 1.0; }
 
@@ -262,6 +266,7 @@ if ($genBamF) {
     $command .= " -b $samInpType $inpF";
     if ($fn_list ne "") { $command .= " 1 $fn_list"; }
     else { $command .= " 0"; }
+    if ($sampling) { $command .= " --sampling"; }
 }
 if ($calcCI) { $command .= " --gibbs-out"; }
 if ($quiet) { $command .= " -q"; }
@@ -478,13 +483,17 @@ Number of threads to use. Both Bowtie and expression estimation will use this ma
 
 Generate a BAM file, 'sample_name.bam', with alignments mapped to genomic coordinates and annotated with their posterior probabilities. In addition, RSEM will call samtools (included in RSEM package) to sort and index the bam file. 'sample_name.sorted.bam' and 'sample_name.sorted.bam.bai' will be generated. (Default: off)
 
+=item B<--sampling-for-bam>
+
+When RSEM generates a BAM file, instead of outputing all alignments a read has with their posterior probabilities, one alignment is sampled and outputed according to the posterior probabilities. If the sampling result is that the read comes from the "noise" transcript, nothing is outputed. It cannot be specified unless --out-bam is specified. (Default: off)
+ 
 =item B<--calc-ci>
 
 Calculate 95% credibility intervals and posterior mean estimates.  (Default: off)
 
 =item B<--seed-length> <int>
 
-Seed length used by the read aligner.  Providing the correct value for this parameter is important for RSEM's accuracy if the data are single-end reads.  If RSEM runs Bowtie, it uses this value for Bowtie's seed length parameter.  (Default: 25)
+Seed length used by the read aligner.  Providing the correct value for this parameter is important for RSEM's accuracy if the data are single-end reads.  If RSEM runs Bowtie, it uses this value for Bowtie's seed length parameter. The minimum value is 25. Any read with its or at least one of its mates' (for paired-end reads) length less than 25 will be ignored. (Default: 25)
 
 =item B<--tag> <string>
 
@@ -624,7 +633,7 @@ Only generated when --out-bam is specified.
 genomic coordinates. Alignments of reads that have identical genomic
 coordinates (i.e., alignments to different isoforms that share the
 same genomic region) are collapsed into one alignment.  The MAPQ field
-of each alignment is set to max(100, floor(-10 * log10(1.0 - w) +
+of each alignment is set to min(100, floor(-10 * log10(1.0 - w) +
 0.5)), where w is the posterior probability of that alignment being
 the true mapping of a read.  In addition, RSEM pads a new tag
 ZW:f:value, where value is a single precision floating number
@@ -694,4 +703,3 @@ Assume the path to the bowtie executables is in the user's PATH environment vari
 
 =cut
 
-#  LocalWords:  usr