Nested Multi-Level Model for RNA-Seq Differential Expression

[jbio]

I have read several other postings regarding nested models, but they did not seem to exactly capture my particular case, and I'm a bit unsure how to proceed with analysis of my model. Any help would be appreciated.

I have a RNA-Seq read count table with (1) multiple tissues from (2) two species and (3) three replicate individuals from each species.

In only three of the tissues from only one species, there is a particular Trait we are interested in. The remainder of the tissues from that species and all the tissues from the other species do not exhibit the 'Trait'.
(Here 'Trait' is a two-level categorical variable). Therefore, the 'Tissue' categories are nested within 'Trait'/non-'Trait'.

The first analysis I tried was just a simple flat comparison with all tissues that exhibit the 'Trait', compared against those that do not exhibit 'Trait.' This, however, falls prey to the classic fallacy when a nested model should be used, and the flat comparison has inflated P-values. This also led to some situations where very high differentials in one Tissue with 'Trait' would compensate for another Tissue with 'Trait' and no difference at all to the non-'Trait' tissues.

So far I have been using limma with voom and have tried several models including:

(1) design <- model.matrix(~0+Trait) using only the 'Trait' coefficient (flat)

(2) design <- model.matrix(~0+Trait/Tissue) using the 'Trait' and all estimable interaction coefficients

(3) design <- model.matrix(~0+Tissue) using all pairwise contrasts between 'Trait and 'nonTrait' to calculate an F-statistic

(3) design <- model.matrix(~0+Tissue) using the contrast (A+B+C)/3 - (D+E+F+G+H+I)/5 where ABC are the 'Trait' and DEFGHI are non-'Trait'

I have also tried these models with/without the duplicateCorrelation() function.

What I want is to have a list where (1) all the tissues with 'Trait' and (2) all the replicates within that tissue both consistently show a differential with all non-'Trait' tissues. The flat comparison is the closest, but does not use the 'Tissue' information to make sure all Tissues have high differentials.

I know from previous postings that most diff. expression software does not explicitly handle mixed or nested models, but if anyone has any tips or has dealt with a similar model, I would appreciate the help. Also, even though I'm currently using limma, I'm perfectly happy to switch if another software will accommodate this model better.

Thanks,

J

[dpryan]

Most DE packages (limma, edgeR, DESeq2) already use shrinkage, which is also the benefit to using a mixed-effect model. Given that, the most obvious model would simply be ~Trait+Tissue, since it sounds like that's what you actually care about.

[jbio]

For topTable then, what coefficients would I want to use?

[dpryan]

Just "Trait".