Comparing Two Conditions Consisting Of Paired Tumor-Normal Samples Using Edger / Limma

Question

[Bontus]

Dear all,

I am wondering how to compare two separate groups that each contain paired tumor normal samples, so assuming we have this experimental setup:

Sample 1_normal      good_outcome
Sample 1_tumor       good_outcome
Sample 2_normal      good_outcome
Sample 2_tumor       good_outcome
Sample 3_normal      bad_outcome
Sample 3_tumor       bad_outcome
Sample 4_normal      bad_outcome
Sample 4_tumor       bad_outcome

is there a way to compare the tumor samples of both groups, taking into account the paired normal samples? My first guess would have been a paired t-test, however, since this is not a before-after comparison I am not sure if it can be applied. To clarify, we have RNA-seq samples with this setup available and I would to use egedR or voom/limma to do this comparison, but I do not know how to create an appropriate contrast matrix (if this is even possible).

Any help is greatly appreciated.

Answer 1

[Devon Ryan]

I'd be surprised if there isn't an example similar to this in either the limma or edgeR user guides or tutorials. In general, you want to apply a model of the form:

count ~ Sample + TumorStatus + Outcome

Have a read through the various user guides, they'll answer most of your questions. I should note that Sample here would be "Sample1","Sample2", etc., since the tumor/normal status has its own factor.

[Chris Cabanski]

See section 3.5 of the edgeR tutorial - this seems to describe the same problem.

[Bontus]

Thank you for your replies. I had a look at the contrast matrix and it appears correctly set up, however if I run the new code as suggested I end up with:

Error in beta[k, ] <- betaj[decr, ] : NAs are not allowed in subscripted assignments

when trying to perform estimateGLMCommonDisp(). Any ideas what might happen here?

[Devon Ryan]

Have a look at this thread on the bioconductor email list: https://stat.ethz.ch/pipermail/bioconductor/2012-April/044811.html (or even this one: https://stat.ethz.ch/pipermail/bioconductor/2012-July/046967.html ). It seems that the model fit just isn't converging, possibly because dispersions are too high. In a few other threads on that email list, people seem to have had success by subsetting their data. Perhaps that will help. If not, you'll likely have to ask Gordon Smyth.