How To Get The Reaction Norm For Each Gene In The Transcriptome?

Question

[johnstantongeddes]

I'm beginning a RNAseq analysis looking for genes that vary in expression across a thermal gradient (12 temperature points) in a non-model (no genome) organism. My data consists of Illumina 100bp PE data for 12 temperature points. Basically, what I'd like to get is the reaction norm for each gene in the transcriptome. This seems like a straight-forward approach to me, but my searching of the literature has not turned up any studies that have done this(!?). To my understanding, the standard methods for gene expression analysis (DeSEQ) only work for pair-wise or multi-factor designs, not a regression design as I've proposed. The closest I've found is this article Abundant gene-by-environment interactions in gene expression reaction norms to copper within Sacharomyces cerevisiae, but this paper uses spotted microarrays (I have Illumina) and thus performs pair-wise comparisons.

I have plans on how to go about this analysis, but am interested in suggestions from the Biostars community or recommendations of citations I may have missed.

Thanks for any suggestions in advance, John

Answer 1

[Steve Lianoglou]

First, could you please add some explaination to your question to clarify what a "reaction norm" is in this context?

In any case, there are ways to analyze time course expression data in a GLM framework using splines. The limma user's guide has an example of how you might do so using microarray data. You might be most comfortable simply shooting your count data through the limma::voom function and proceed as prescribed there.

[Neilfws]

I assume a reaction norm is the same thing as a norm of reaction :)

[johnstantongeddes]

Yes, thank you Neilfws, I am referring to a norm of reaction. I'm looking over the limma users's guide...the "time series analysis" is the type of analysis I was looking for, so I'll likely accept this answer but am curious about other approaches out there...

[Steve Lianoglou]

I actually stumbled on that page when I googled for "reaction norm" before posting this answer. I still don't understand how you'd use "that" in this context. Do you just want to estimate a gaussian per gene that best fits its expression across all 12 temp time points? If so, then what?

[johnstantongeddes]

I have two samples that I a priori expect to differ in thermal tolerance. The idea is to fit a polynomial to expression across the 12 temp points for each sample. Obviously you can always fit a polynomial, so I'll randomize the expression values among treatments many times and only consider significant the genes with linear and quadratic terms in the tails of the empirical distributions. Then identify genes that differ in intercept, linear and quadratic terms among the two samples. I'm still working through the limma guide/examples so not sure if it will work with this, or I'll just code it myself.