<p>I am looking for input on ways to put numbers on how much different factors/covariates contribute to gene expression patterns. For example, let's say we have RNA-seq (or microarray) data from male and female mice that have received two different treatments (or no treatment) and the samples have been prepped in two different batches. Now we would like a bar plot (or similar) of how much the different covariates (1) sex, (2) treatment and (3) batch effects contribute, respectively, to the gene expression. </p>

<p>Some ideas:</p>

<ol>
<li>Use PCA on the expression matrix and calculate the correlation between PC1 scores and the covariate vectors.</li>
<li>Use some machine learning algorithm to try to discriminate between male/female (or treatments, or batches), and evaluate by cross-validation how well this is working. Use the prediction performance as a measure of how strongly the covariate is reflected in the expression profile.</li>
<li>Is there a way to find out this sort of summarized importance from the linear models fitted as part of limma/edgeR/DESeq workflows? </li>
</ol>

<p>I'd be grateful for any thoughts.</p>

<p>As a first step, it might be useful to remove genes with count 0 in more than 90% of your samples. Then obtain normalized and log-transformed expression estimates (logFPKM; maybe edgeR and voom() can be of use?) for all your samples and genes. I assume you already have obtained them somehow. Then prepare your data in such a way that you have a data frame where a row represents 1 measurement and the columns geneId, sex, treatment, batch, logFpkm (maybe aggregate() or melt() functions can be useful to reshape your data). Then perform anova and look at the sums of squares and the p-values of the F-test. Most likely, geneId has the biggest influence on your logFpkms and hopefully treatment is a good second.</p>

<pre><code>fit &lt;- lm(logFpkm ~  geneId + sex + treatment + batch,data=yourData)
anova(fit)
</code></pre>

<p>When you have many samples and many genes you can also take a random sample of your genes (say 1000) to get an idea and then scale up to bigger sample sizes.</p>

Edit: Brentp, you are right. Its maybe better to exclude the geneId from the linear model
