7
$\begingroup$

Without going into too much background, I just joined up with a lab as a bioinformatics intern while I'm completing my masters degree in the field. The lab has data from an RNA-seq they outsourced, but the only problem is that the only data they have is preprocessed from the company that did the sequencing: filtering the reads, aligning them, and putting the aligned reads through RSEM. I currently have output from RSEM for each of the four samples consisting of: gene id, transcript id(s), length, expected count, and FPKM. I am attempting to get the FASTQ files from the sequencing, but for now, this is what I have, and I'm trying to get something out of it if possible.

I found this article that talks about how expected read counts can be better than raw read counts when analyzing differential expression using EBSeq; it's just one guy's opinion, and it's from 2014, so it may be wrong or outdated, but I thought I'd give it a try since I have the expected counts.

However, I have just a couple of questions about running EBSeq that I can't find the answers to:

1: In the output RSEM files I have, not all genes are represented in each, about 80% of them are, but for the ones that aren't, should I remove them before analysis with EBSeq? It runs when I do, but I'm not sure if it is correct.

2: How do I know which normalization factor to use when running EBSeq? This is more of a conceptual question rather than a technical question.

Thanks!

Improve this question
$\endgroup$

2 Answers 2

Reset to default
6
$\begingroup$

Yes, that blog post does represent just one guy's opinion (hi!) and it does date all the way back to 2014, which is, like, decades in genomics years. :-) By the way, there is quite a bit of literature discussing the improvements that expected read counts derived from an Expectation Maximization algorithm provide over raw read counts. I'd suggest reading the RSEM papers for a start[1][2].

But your main question is about the mechanics of running RSEM and EBSeq. First, RSEM was written explicitly to be compatible with EBSeq, so I'd be very surprised if it does not work correctly out-of-the-box. Second, EBSeq's MedianNorm function worked very well in my experience for normalizing the library counts. Along those lines, the blog you mentioned above has another post that you may find useful.

But all joking aside, these tools are indeed dated. Alignment-free RNA-Seq tools provide orders-of-magnitude improvements in runtime over the older alignment-based alternatives, with comparable accuracy. Sailfish was the first in a growing list of tools that now includes Salmon and Kallisto. When starting a new analysis from scratch (i.e. if you ever get the original FASTQ files), there's really no good reason not to estimate expression using these much faster tools, followed by a differential expression analysis with DESeq2, edgeR, or sleuth.

1Li B, Ruotti V, Stewart RM, Thomson JA, Dewey CN (2010) RNA-Seq gene expression estimation with read mapping uncertainty. Bioinformatics, 26(4):493–500, doi:10.1093/bioinformatics/btp692.

2Li B, Dewey C (2011) RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics, 12:323, doi:10.1186/1471-2105-12-323.

Improve this answer
$\endgroup$
5
  • 2
    $\begingroup$ "There's really no good reason not to estimate expression using these much faster tools" -- unless you don't have the raw reads, as is the case here $\endgroup$
    – gringer ♦
    Jun 2, 2017 at 1:32
  • $\begingroup$ Oh wow. Big oversight on my part! $\endgroup$
    – Daniel Standage
    Jun 2, 2017 at 1:58
  • $\begingroup$ Wow I never expected to get a response from the actual author! I did use R to process the data frames into a single matrix of all the expected counts for each gene of each sample. $\endgroup$
    – J0HN_TIT0R
    Jun 2, 2017 at 17:04
  • $\begingroup$ I was just as surprised to see a link to my old blog on StackExchange! :) $\endgroup$
    – Daniel Standage
    Jun 2, 2017 at 17:07
  • $\begingroup$ As a small follow-up question, I'm ultimately trying to obtain the fold change and associated p-value for each gene per condition. I've found GetMultiFC() to get the fold changes, but I'm unclear about the exact difference between the fold change and the posterior fold change. I'm thinking that the posterior fold change is just the fold change for the normalized values, so it's the one I should be using, but I'm not sure about that. Also, is it possible to get associated p-values? $\endgroup$
    – J0HN_TIT0R
    Jun 2, 2017 at 19:26
3
$\begingroup$
  1. Include all genes/transcripts in your analysis.

A transcript that is not detected could be undetected through sampling error (i.e. the sequencer / library prep just happened to miss that transcript), or it could be because the transcript isn't generated in a particular sample. It's not uncommon for genes to be switched off in response to different biological factors, so zero-count genes shouldn't be ignored. I can't speak from experience with EBSeq, but as long as the analysis package treats a zero count as "unobserved" rather than "absent" (and makes relevant corrections), it's a good idea to keep them in.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.