#const
my $BURNIN = 200;
-my $CHAINLEN = 1000;
+my $NCV = 1000;
my $SAMPLEGAP = 1;
my $CONFIDENCE = 0.95;
my $NSPC = 50;
my $genGenomeBamF = 0;
my $sampling = 0;
my $calcCI = 0;
+my $var_opt = 0; # temporarily, only for internal use
my $quiet = 0;
my $help = 0;
"estimate-rspd" => \$estRSPD,
"num-rspd-bins=i" => \$B,
"p|num-threads=i" => \$nThreads,
+ "no-bam-output" => sub { $genBamF = 0; },
"output-genome-bam" => \$genGenomeBamF,
"sampling-for-bam" => \$sampling,
+ "var" => \$var_opt,
"calc-ci" => \$calcCI,
"ci-memory=i" => \$NMB,
"time" => \$mTime,
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 5!\n", -exitval => 2, -verbose => 2) if ($L < 5);
-pod2usage(-msg => "--sampling-for-bam cannot be specified if --out-bam is not specified!\n", -exitval => 2, -verbose => 2) if ($sampling && !$genBamF);
+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);
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"; }
$imdName = "$temp_dir/$sampleToken";
-if (!$is_sam && !$is_bam && $phred33 + $phred64 + $solexa == 0) { $phred33 = 1; }
+if (!$is_sam && !$is_bam && !$no_qual && ($phred33 + $phred64 + $solexa == 0)) { $phred33 = 1; }
my ($mate_minL, $mate_maxL) = (1, $maxL);
if (!$is_sam && !$is_bam) {
$command = $bowtie_path."bowtie";
- if ($read_type == 0 || $read_type == 2) { $command .= " -f"; }
+ if ($no_qual) { $command .= " -f"; }
else { $command .= " -q"; }
if ($phred33) { $command .= " --phred33-quals"; }
elsif ($phred64) { $command .= " --phred64-quals"; }
elsif ($solexa) { $command .= " --solexa-quals"; }
- else { print "Oh, no!!!"; exit(2); }
$command .= " -n $C -e $E -l $L";
if ($read_type == 2 || $read_type == 3) { $command .= " -I $minL -X $maxL"; }
$command .= " -1 $mate1_list -2 $mate2_list";
}
- $command .= " | gzip > $imdName.sam.gz";
+ $command .= " | gzip > $sampleName.sam.gz";
if ($mTime) { $time_start = time(); }
if ($mTime) { $time_end = time(); $time_alignment = $time_end - $time_start; }
- $inpF = "$imdName.sam.gz";
+ $inpF = "$sampleName.sam.gz";
$is_sam = 1; # output of bowtie is a sam file
}
if ($mTime) { $time_start = time(); }
-if ($calcCI) {
- $command = $dir."rsem-run-gibbs $refName $sampleName $sampleToken $BURNIN $CHAINLEN $SAMPLEGAP";
-# $command .= " -p $nThreads";
+if ($calcCI || $var_opt) {
+ $command = $dir."rsem-run-gibbs $refName $sampleName $sampleToken $BURNIN $NCV $SAMPLEGAP";
+ $command .= " -p $nThreads";
+ if ($var_opt) { $command .= " --var"; }
if ($quiet) { $command .= " -q"; }
&runCommand($command);
+}
+if ($calcCI) {
system("mv $sampleName.isoforms.results $imdName.isoforms.results.bak1");
system("mv $sampleName.genes.results $imdName.genes.results.bak1");
&collectResults("$imdName.iso_res", "$sampleName.isoforms.results"); # isoform level
&collectResults("$imdName.gene_res", "$sampleName.genes.results"); # gene level
- $command = $dir."rsem-calculate-credibility-intervals $refName $sampleName $sampleToken $CONFIDENCE $NSPC $NMB";
+ $command = $dir."rsem-calculate-credibility-intervals $refName $sampleName $sampleToken $CONFIDENCE $NCV $NSPC $NMB";
+ $command .= " -p $nThreads";
if ($quiet) { $command .= " -q"; }
&runCommand($command);
if ($mTime) {
open(OUTPUT, ">$sampleName.time");
- print OUTPUT "Alignment: $time_alignment s.\n";
- print OUTPUT "RSEM: $time_rsem s.\n";
- print OUTPUT "CI: $time_ci s.\n";
+ print OUTPUT "Aligning reads: $time_alignment s.\n";
+ print OUTPUT "Estimating expression levels: $time_rsem s.\n";
+ print OUTPUT "Calculating credibility intervals: $time_ci s.\n";
my $time_del = $time_end - $time_start;
- print OUTPUT "Delete: $time_del s.\n";
+# print OUTPUT "Delete: $time_del s.\n";
close(OUTPUT);
}
Number of threads to use. Both Bowtie and expression estimation will use this many threads. (Default: 1)
+=item B<--no-bam-output>
+
+Do not output any BAM file. (Default: off)
+
=item B<--output-genome-bam>
Generate a BAM file, 'sample_name.genome.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.genome.sorted.bam' and 'sample_name.genome.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. (Default: off)
-
+
=item B<--calc-ci>
Calculate 95% credibility intervals and posterior mean estimates. (Default: off)
=item B<--ci-memory> <int>
-Amount of memory (in MB) RSEM is allowed to use for computing credibility intervals. (Default: 1024)
+Maximum size (in memory, MB) of the auxiliary buffer used for computing credibility intervals (CI). Set it larger for a faster CI calculation. However, leaving 2 GB memory free for other usage is recommended. (Default: 1024)
=item B<--keep-intermediate-files>
Keep temporary files generated by RSEM. RSEM creates a temporary directory, 'sample_name.temp', into which it puts all intermediate output files. If this directory already exists, RSEM overwrites all files generated by previous RSEM runs inside of it. By default, after RSEM finishes, the temporary directory is deleted. Set this option to prevent the deletion of this directory and the intermediate files inside of it. (Default: off)
+=item B<--time>
+
+Output time consumed by each step of RSEM to 'sample_name.time'. (Default: off)
+
=item B<-q/--quiet>
Suppress the output of logging information. (Default: off)
=item B<sample_name.transcript.bam, sample_name.transcript.sorted.bam and sample_name.transcript.sorted.bam.bai>
+Only generated when --no-bam-output is not specified.
+
'sample_name.transcript.bam' is a BAM-formatted file of read
alignments in transcript coordinates. The MAPQ field of each alignment
is set to min(100, floor(-10 * log10(1.0 - w) + 0.5)), where w is the
=item B<sample_name.genome.bam, sample_name.genome.sorted.bam and sample_name.genome.sorted.bam.bai>
-Only generated when --output-genome-bam is specified.
+Only generated when --no-bam-output is not specified and --output-genome-bam is specified.
'sample_name.genome.bam' is a BAM-formatted file of read alignments in
genomic coordinates. Alignments of reads that have identical genomic
'sample_name.genome.sorted.bam' and 'sample_name.genome.sorted.bam.bai' are the
sorted BAM file and indices generated by samtools (included in RSEM package).
+=item B<sample_name.sam.gz>
+
+Only generated when the input files are raw reads instead of SAM/BAM format files
+
+It is the gzipped SAM output produced by bowtie aligner.
+
+=item B<sample_name.time>
+
+Only generated when --time is specified.
+
+It contains time (in seconds) consumed by aligning reads, estimating expression levels and calculating credibility intervals.
+
=item B<sample_name.stat>
This is a folder instead of a file. All model related statistics are stored in this folder. Use 'rsem-plot-model' can generate plots using this folder.
projects / rsem.git / blobdiff