my $genGenomeBamF = 0;
my $sampling = 0;
my $calcCI = 0;
+my $var_opt = 0; # temporarily, only for internal use
my $quiet = 0;
my $help = 0;
my $mTime = 0;
my ($time_start, $time_end, $time_alignment, $time_rsem, $time_ci) = (0, 0, 0, 0, 0);
+my $mate1_list = "";
+my $mate2_list = "";
+my $inpF = "";
+
+my ($refName, $sampleName, $sampleToken, $temp_dir, $stat_dir, $imdName, $statName) = ();
+my $gap = 32;
+
GetOptions("keep-intermediate-files" => \$keep_intermediate_files,
+ "temporary-folder=s" => \$temp_dir,
"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,
+ "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"; }
if ($strand_specific) { $probF = 1.0; }
-my $mate1_list = "";
-my $mate2_list = "";
-my $inpF = "";
-
-my ($refName, $sampleName, $sampleToken, $temp_dir, $stat_dir, $imdName) = ();
-my $gap = 32;
-
if ($paired_end) {
if ($no_qual) { $read_type = 2; }
else { $read_type = 3; }
if ($pos < 0) { $sampleToken = $sampleName; }
else { $sampleToken = substr($sampleName, $pos + 1); }
-$temp_dir = "$sampleName.temp";
+if ($temp_dir eq "") { $temp_dir = "$sampleName.temp"; }
$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";
+$statName = "$stat_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 > $sampleName.sam.gz";
+ # pipe to samtools to generate a BAM file
+ $command .= " | $dir\sam/samtools view -S -b -o $imdName.bam -";
if ($mTime) { $time_start = time(); }
if ($mTime) { $time_end = time(); $time_alignment = $time_end - $time_start; }
- $inpF = "$sampleName.sam.gz";
- $is_sam = 1; # output of bowtie is a sam file
+ $inpF = "$imdName.bam";
+ $is_bam = 1; # alignments are outputed as a BAM file
}
if ($mTime) { $time_start = time(); }
-$command = $dir."rsem-parse-alignments $refName $sampleName $sampleToken";
+$command = $dir."rsem-parse-alignments $refName $imdName $statName";
my $samInpType;
if ($is_sam) { $samInpType = "s"; }
print OUTPUT "$L\n";
close(OUTPUT);
-$command = $dir."rsem-run-em $refName $read_type $sampleName $sampleToken -p $nThreads";
+$command = $dir."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 ($calcCI) { $command .= " --gibbs-out"; }
+if ($calcCI || $var_opt) { $command .= " --gibbs-out"; }
if ($quiet) { $command .= " -q"; }
&runCommand($command);
+&collectResults("$imdName.iso_res", "$sampleName.isoforms.results"); # isoform level
+&collectResults("$imdName.gene_res", "$sampleName.genes.results"); # gene level
+
if ($genBamF) {
$command = $dir."sam/samtools sort $sampleName.transcript.bam $sampleName.transcript.sorted";
&runCommand($command);
}
}
-&collectResults("$imdName.iso_res", "$sampleName.isoforms.results"); # isoform level
-&collectResults("$imdName.gene_res", "$sampleName.genes.results"); # gene level
-
if ($mTime) { $time_end = time(); $time_rsem = $time_end - $time_start; }
if ($mTime) { $time_start = time(); }
-if ($calcCI) {
- $command = $dir."rsem-run-gibbs $refName $sampleName $sampleToken $BURNIN $NCV $SAMPLEGAP";
+if ($calcCI || $var_opt) {
+ $command = $dir."rsem-run-gibbs $refName $imdName $statName $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
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)
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<--temporary-folder> <string>
+
+Set where to put the temporary files generated by RSEM. If the folder specified does not exist, RSEM will try to create it. (Default: sample_name.temp)
+
=item B<--time>
Output time consumed by each step of RSEM to 'sample_name.time'. (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
projects / rsem.git / blobdiff