I am investigating the effect of statins on the expression of approximately 50 000 genes by comparing two groups: patients treated with or without statins.
I chose a model as follows: gene expression ~ age + gender + smoking status. Unfortunately, after adjustment, the qqplot of p-values obtained after comparing the two groups (Wilcoxon signed rank test) indicates a bias. I have found methods that could help to improve the model (AIC, Mallows Cp, R square), however they were based on the expression of a single gene.
The method "surrogate variable analysis" (SVA) has the aim to remove factors that are unknown, unmeasured or not taken in account in a model and takes in account all gene expressions. However, I have some difficulties to make it work on my data.
Do you know a method that allows to test a model on a set of genes and other methods similar to SVA ?
Hi Jerome. Good question. Can you tell us more about what software you are using for the analysis? Maybe it could help people to provide you with a more usable answer.
R is the programming language that I used to test the methods (package for rlm, AIC and sva). Note that there is a web application of the sva method: http://psychiatry.igm.jhmi.edu/sva/, I could not use it because the data to analyse is too big.
Providing the code you used would be helpful, but also: are you sure you need to use sva? Do you expect strong batch effects in the data? If you do a PCA plot of your arrays, do they cluster in strange ways?
I managed to make the standalone version of SVA work on my data. The qqplot looks better but there is still a bias. I use SVA because I do not know which variable(s) is (are) causing the bias and if the bias has been measured in my variable(s). I do not know what to expect about the batch effect in the data, so I guess in this case the best thing to do is to act as if there were some and SVA seems be suited for this. I am working on the PCA plot to see if there is anything strange.
There is something I do not understand in your model. You said you want to test the effect of statins on the expression of genes. However statins are not included in your model. How can you then test for statins effect? Or is this a different model used in your project for statins?