my $mTime = 0;
my ($time_start, $time_end, $time_alignment, $time_rsem, $time_ci) = (0, 0, 0, 0, 0);
my $mTime = 0;
my ($time_start, $time_end, $time_alignment, $time_rsem, $time_ci) = (0, 0, 0, 0, 0);
"no-qualities" => \$no_qual,
"paired-end" => \$paired_end,
"strand-specific" => \$strand_specific,
"no-qualities" => \$no_qual,
"paired-end" => \$paired_end,
"strand-specific" => \$strand_specific,
"estimate-rspd" => \$estRSPD,
"num-rspd-bins=i" => \$B,
"p|num-threads=i" => \$nThreads,
"estimate-rspd" => \$estRSPD,
"num-rspd-bins=i" => \$B,
"p|num-threads=i" => \$nThreads,
"output-genome-bam" => \$genGenomeBamF,
"sampling-for-bam" => \$sampling,
"output-genome-bam" => \$genGenomeBamF,
"sampling-for-bam" => \$sampling,
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 => "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"; }
if ($strand_specific) { $probF = 1.0; }
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"; }
if ($strand_specific) { $probF = 1.0; }
if ($pos < 0) { $sampleToken = $sampleName; }
else { $sampleToken = substr($sampleName, $pos + 1); }
if ($pos < 0) { $sampleToken = $sampleName; }
else { $sampleToken = substr($sampleName, $pos + 1); }
$stat_dir = "$sampleName.stat";
if (!(-d $temp_dir) && !mkdir($temp_dir)) { print "Fail to create folder $temp_dir.\n"; exit(-1); }
if (!(-d $stat_dir) && !mkdir($stat_dir)) { print "Fail to create folder $stat_dir.\n"; exit(-1); }
$imdName = "$temp_dir/$sampleToken";
$stat_dir = "$sampleName.stat";
if (!(-d $temp_dir) && !mkdir($temp_dir)) { print "Fail to create folder $temp_dir.\n"; exit(-1); }
if (!(-d $stat_dir) && !mkdir($stat_dir)) { print "Fail to create folder $stat_dir.\n"; exit(-1); }
$imdName = "$temp_dir/$sampleToken";
else { $command .= " -q"; }
if ($phred33) { $command .= " --phred33-quals"; }
elsif ($phred64) { $command .= " --phred64-quals"; }
elsif ($solexa) { $command .= " --solexa-quals"; }
else { $command .= " -q"; }
if ($phred33) { $command .= " --phred33-quals"; }
elsif ($phred64) { $command .= " --phred64-quals"; }
elsif ($solexa) { $command .= " --solexa-quals"; }
$command .= " -n $C -e $E -l $L";
if ($read_type == 2 || $read_type == 3) { $command .= " -I $minL -X $maxL"; }
$command .= " -n $C -e $E -l $L";
if ($read_type == 2 || $read_type == 3) { $command .= " -I $minL -X $maxL"; }
if ($genBamF) {
$command .= " -b $samInpType $inpF";
if ($fn_list ne "") { $command .= " 1 $fn_list"; }
else { $command .= " 0"; }
if ($sampling) { $command .= " --sampling"; }
}
if ($genBamF) {
$command .= " -b $samInpType $inpF";
if ($fn_list ne "") { $command .= " 1 $fn_list"; }
else { $command .= " 0"; }
if ($sampling) { $command .= " --sampling"; }
}
if ($mTime) { $time_end = time(); $time_rsem = $time_end - $time_start; }
if ($mTime) { $time_start = time(); }
if ($mTime) { $time_end = time(); $time_rsem = $time_end - $time_start; }
if ($mTime) { $time_start = time(); }
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
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
- if ($status != 0) {
- my $errmsg;
- if (scalar(@_) > 1) { $errmsg = $_[1]; }
- else { $errmsg = "\"$command\" failed! Plase check if you provide correct parameters/options for the pipeline!"; }
- print $errmsg."\n";
- exit(-1);
+ if ($status != 0) {
+ my $errmsg = "";
+ if (scalar(@_) > 1) { $errmsg .= $_[1]."\n"; }
+ $errmsg .= "\"$_[0]\" failed! Plase check if you provide correct parameters/options for the pipeline!\n";
+ print $errmsg;
+ exit(-1);
=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>
=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)
+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<--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<--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)
'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
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 representing the posterior
'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
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 representing the posterior
+probability. Because this file contains all alignment lines produced
+by bowtie or user-specified aligners, it can also be used as a
+replacement of the aligner generated BAM/SAM file. For paired-end
+reads, if one mate has alignments but the other does not, this file
+marks the alignable mate as "unmappable" (flag bit 0x4) and appends an
+optional field "Z0:A:!".
'sample_name.transcript.sorted.bam' and
'sample_name.transcript.sorted.bam.bai' are the sorted BAM file and
'sample_name.transcript.sorted.bam' and
'sample_name.transcript.sorted.bam.bai' are the sorted BAM file and
'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.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).
'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).