Question on expected counts/wiggle plots

[Simon Moxon]

Hi,

I'm new to using RSEM and am a little confused by the results I'm getting (I'm using version 1.2.4). I've followed the instructions on the website and given a transcript fasta file to rsem-prepare-reference to create the necessary files. Next I ran rsem-calculate-expression on several samples using the prepared transcript file with paired-end reads - this all works fine and I get the output files. I then used the *.isoforms.results to calculate differential expression using DESeq (also fine). However I then used rsem-plot-transcript-wiggles to generate pdfs for my differentially expressed transcripts and I was a bit surprised to see that they differ wildly from the bam file output (generated using rsem-calculate-expression). I've attached screenshots of the plots of a particular transcript in two samples along with a genome browser plot of the same gene (using the bam (blue) & wig (black) files generated by RSEM).

The output from rsem-calculate-expression for the transcript in question is shown below - counts are highlighted:

NC_untreated:

gi|148234886|ref|NM_001090653| gi|148234886|ref|NM_001090653| 2487 2375.77 338.30 8.81 6.50 100.00

NC_treated:

gi|148234886|ref|NM_001090653| gi|148234886|ref|NM_001090653| 2487 2363.61 0.08 0.00 0.00 100.00

As you can see the coverage on the bam tracks (genome browser) are similar but the expected counts are very different for this transcript in the two samples. Can anyone explain to me why this is?

Thanks

Simon

NC_untreated

NC_treated

[Colin Dewey]

Hi Simon,

Noting the differences in the y-axis scales (NC_untreated has a max of 200 or so, whereas NC_treated has a max of ~0.04), it looks like the plots agree with the counts.  Does that address your question?

Best,

Colin

--
RSEM website: http://deweylab.biostat.wisc.edu/rsem/
---
You received this message because you are subscribed to the Google Groups "RSEM Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to rsem-users+...@googlegroups.com.
To post to this group, send email to rsem-...@googlegroups.com.
Visit this group at http://groups.google.com/group/rsem-users?hl=en-US.
 
 

[Simon Moxon]

Hi Colin,

Thanks for your reply. Yes the plots definitely agree with the counts. I'm just not sure why the counts for NC_treated were so low. Looking at the bam file tracks the coverage looks similar across the transcript in the two samples (albeit slightly lower in treated vs untreated) however the RSEM counts are significantly different. That's why I was a bit surprised.

Thanks

Simon

[Colin Dewey]

Hi Simon,

Sorry, I had not noticed your genome browser plot.  First, a couple of things to note about the bam and wig files generated by RSEM (1) the bam file contains all candidate alignments for each read with each alignment annotated with a posterior probability of it being the "true" alignment (2) the wig files are generated by summing the *posterior probabilities* of the reads overlapping each genomic position and not by simply counting the number of reads with alignments at that position.  My guess is that your genome browser is generating wiggle plots from the BAM by simply counting the number of alignments at each position, without regard to their posterior probabilities.  In addition, I'm guessing that this gene has a close paralog somewhere else in the genome. In the untreated sample RSEM is assigning most of the reads to this gene, whereas in the treated sample RSEM is assigning most of the reads to the paralog.

Colin

[Simon Moxon]

Thanks for the explanation Colin. I was guessing that it was due the bulk of the reads present in the bam alignment being assigned to another transcript but just wanted to check that this was the most likely explanation.

Thanks again for your help with this - it's much appreciated!

Simon