Hi All, 

Thanks in advance for the help! I am still new to edgeR and welcome for any suggestions!

I am trying to figure out which method is better off to analyze my data. I am interesting in translational change before vs. after stress in yeast(s. cerevisiae). I carried out polysome profiles before and after stress, then RNA-seq of collected fractions of whole polysome profiles (total 9 fractions per polysome profile). During RNA sample extraction, I spike-in equal amount of RNA (from S.pombe) for later normalization. I used the normalized counts (based on spike-in counts), summed up the reads of polysome fractions(poly) and reads of the whole polysome profile(total), and tried to use the ratio (poly/total) to access the translation change. I have 2 biological replicates and 2 conditions, before and after stress. I have used the ratio number in limma and got reasonable results. Meanwhile, I also get suggestions to use edgeR with summed poly-counts as input and total counts as offsets. But, when I run the edgeR, it crushed after a certain code. Here is my code for running edgeR:

> x = DGEList(counts = inp) # inp contains normalized reads of summed polysome fractions
> x$offest = off # off is the sum of total reads, have the same dimension as inp
> x = estimateCommonDisp(x)

>x = estimateGLMTrendedDisp(x, design)

The R is crushed on this step.

Any suggestions and thoughts with my codes? I am also wondering what people think about using limma vs. edgeR for this analysis.

Thanks,

Elisha

To summarize your experimental design and data:

• You have 2 conditions, before and after stress.
• Each condition has two biological replicates.
• Each replicate has two counts - ribosome-bound and total - for each gene.

I would not use the total counts to compute offsets here, as edgeR treats the supplied offsets as known values. This means that variability in the total counts is not properly modelled. Instead, I would do something like this:

condition <- rep(c("before", "after"), each=4)
replicate <- gl(4, 2)
count <- rep(c("ribo", "total"), 4)
design <- model.matrix(~0 + replicate + count:condition)
design <- design[,!grepl("counttotal", colnames(design))] # get to full rank

The first four columns of design represent the log-total expression in each replicate, and are not useful. The second last column (countribo:conditionafter) represents the log-fold change of the ribosomal count over the total count in the after-stress condition, while the last column (countribo:conditionbefore) represents the log-fold change of the ribosomal count over the total count in the before stress condition. It is then simple to test for changes in ribosomal/total ratios between conditions:

colnames(design) <- make.names(colnames(design)) # syntactically valid.
con <- makeContrasts(countribo.conditionafter - countribo.conditionbefore, levels=design)

Edit: The advice above assumes that the total and ribosome-bound counts came from two separate experiments, much like an IP and input control in ChIP-seq experiments. However, if the total count was obtained by adding the ribosome-bound count to a non-ribosome-bound count, then the total and ribosome-bound counts would be correlated. This would not be good as it would invalidate independence assumptions in the model. Instead, you should model the non-ribosome-bound count directly:

count <- rep(c("ribo", "nonribo"), 4)

You may also be interested in the modelMatrixMeth function, which provides a convenient way of setting up the design matrix above. It was initially designed for methylation experiments (methylated vs. unmethylated counts) but is equally applicable here.

Following up on Aaron's answer, you can analyse your bound vs non-bound counts in edgeR as one would analyse methylated counts, see:

http://www.statsci.org/smyth/pubs/edgeRMethylationPreprintJan2018.pdf

Your data is exactly analogous, with your bound counts corresponding to methylated counts and non-bound counts corresponding to unmethylated counts.