Made the BAM output of RSEM have all information contained in the input alignments...

[rsem.git] / rsem-calculate-expression

diff --git a/rsem-calculate-expression b/rsem-calculate-expression
index 074d895ecf173cad4d5b0866c42b463655365747..457f1771e5314c2d29b46a3ab16a2a9e3622e49e 100755 (executable)
--- a/rsem-calculate-expression
+++ b/rsem-calculate-expression
@@ -62,7 +62,15 @@ my $strand_specific = 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,
@@ -126,13 +134,6 @@ if ($L < 25) { print "Warning: the seed length set is less than 25! This is only
 
 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; }
@@ -167,13 +168,14 @@ my $pos = rindex($sampleName, '/');
 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 && !$no_qual && ($phred33 + $phred64 + $solexa == 0)) { $phred33 = 1; }
 
@@ -211,7 +213,8 @@ if (!$is_sam && !$is_bam) {
        $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(); }
 
@@ -219,13 +222,13 @@ if (!$is_sam && !$is_bam) {
 
     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"; } 
@@ -258,7 +261,7 @@ print OUTPUT "$mean $sd\n";
 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"; }
@@ -270,6 +273,9 @@ 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);
@@ -286,15 +292,12 @@ if ($genBamF) {
     }
 }
 
-&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 || $var_opt) {
-    $command = $dir."rsem-run-gibbs $refName $sampleName $sampleToken $BURNIN $NCV $SAMPLEGAP";
+    $command = $dir."rsem-run-gibbs $refName $imdName $statName $BURNIN $NCV $SAMPLEGAP";
     $command .= " -p $nThreads";
     if ($var_opt) { $command .= " --var"; }
     if ($quiet) { $command .= " -q"; }
@@ -561,6 +564,10 @@ Maximum size (in memory, MB) of the auxiliary buffer used for computing credibil
 
 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)