I am attempting to test some hypotheses about whether thermal effects on transcription (polyA RNAseq data) are tissue-specific. The experiment includes four replicate samples per tissue/temperature combination, with three tissues and four temperatures (the treatment).

Ordinarily, I'd test for factor interactions in a model implemented in edgeR. But, this design involves multiple tissues with widespread differences in expression, which I believe does not play well with the standard normalization methods - in edgeR, the calculation of the normalization scaling factors.

The qsmooth normalization method (Hicks, Stephanie C., et al. "Smooth quantile normalization." Biostatistics 19.2 (2018): 185-198.) seems suitable here, but I haven't been able to track down whether data normalized by qsmooth would be appropriate for the edgeR (or deseq2) models. As far as I can tell, most (maybe all?) pre-normalizations screw up the mean-variance relationships assumed in the models.

So, three questions:

1) is there any way to make qsmooth work with edgeR or deseq2 (or voom)? 2) If not, are there any other end-arounds to analyze data including large tissue effects in these packages? 3) Are there other models/packages that would be more appropriate?

Thanks!

edgeR can handle pretty much any normalization method, without messing up the mean-variance relationship, through the use of an offset matrix. I don't know how that would be done with qsmooth though. The easiest way to apply qsmooth I guess would be to use the limma-trend pipeline and apply qsmooth to the logCPM values.

I'd also like to see some testing of qsmooth to show that qsmooth-normalized data still controls the FDR rate correctly. As far as I know, no one has shown that yet.