From 9210a5fece7ec2854eb834d5b2dcbe2d12fbebf1 Mon Sep 17 00:00:00 2001
From: Bo Li <bli@cs.wisc.edu>
Date: Sun, 15 Apr 2012 08:27:02 -0500
Subject: [PATCH] Added some instructions on how to visualize transcript
 coordinate BAM/WIG files using IGV

---
 README.md                 | 12 ++++++++++++
 rsem-calculate-expression |  7 ++++++-
 2 files changed, 18 insertions(+), 1 deletion(-)

diff --git a/README.md b/README.md
index 6e8acd8..f582b33 100644
--- a/README.md
+++ b/README.md
@@ -169,6 +169,18 @@ For UCSC genome browser, please refer to the [UCSC custom track help page](http:
 
 For integrative genomics viewer, please refer to the [IGV home page](http://www.broadinstitute.org/software/igv/home). Note: Although IGV can generate read depth plot from the BAM file given, it cannot recognize "ZW" tag RSEM puts. Therefore IGV counts each alignment as weight 1 instead of the expected weight for the plot it generates. So we recommend to use the wiggle file generated by RSEM for read depth visualization.
 
+Here are some guidance for visualizing transcript coordinate files:
+
+1) Import the transcript sequences as a genome 
+
+Select File -> Import Genome, then fill in ID, Name and Fasta file. Fasta file should be 'reference_name.transcripts.fa'. After that, click Save button. Suppose ID is filled as 'reference_name', a file called 'reference_name.genome' will be generated. Next time, we can use: File -> Load Genome, then select 'reference_name.genome'.
+
+2) Load visualization files
+
+Select File -> Load from File, then choose one transcript coordinate visualization file generated by RSEM. IGV might require you to convert wiggle file to tdf file. You should use igvtools to perform this task. One way to perform the conversion is to use the following command
+
+    igvtools tile reference_name.transcript.wig reference_name.transcript.tdf reference_name.genome   
+ 
 #### c) Generating Transcript Wiggle Plots
 
 To generate transcript wiggle plots, you should run the
diff --git a/rsem-calculate-expression b/rsem-calculate-expression
index c01c70a..2a09c15 100755
--- a/rsem-calculate-expression
+++ b/rsem-calculate-expression
@@ -49,6 +49,7 @@ my $genBamF = 1;  # default is generating transcript bam file
 my $genGenomeBamF = 0;
 my $sampling = 0;
 my $calcCI = 0;
+my $var_opt = 0; # temporarily, only for internal use
 my $quiet = 0;
 my $help = 0;
 
@@ -89,6 +90,7 @@ GetOptions("keep-intermediate-files" => \$keep_intermediate_files,
 	   "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,
@@ -291,12 +293,15 @@ if ($mTime) { $time_end = time(); $time_rsem = $time_end - $time_start; }
 
 if ($mTime) { $time_start = time(); }
 
-if ($calcCI) {
+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
-- 
2.39.2