I have a RNAseq dataset for drug response. The data is voom normalized. The dataset consists of two species. Under each species there are three groups - control, susceptible and resilient. And there are three replicates for each group. Thus the total number of samples is 18.
The samples are represented with following labels:
Species1 (control, susceptible, resilient) <- c("a","a","a","b","b","b","c","c","c")
Species2 (control, susceptible, resilient) <- c("d","d","d","e","e","e","f","f","f")
In order to perform Gene Set Analysis within species (for example, within species1), I do the following:
y <- matrix(rnorm(1000*18),1000,18)
group = as.factor(c("a","a","a","b","b","b","c","c","c","d","d","d","e","e","e","f","f","f"))
des <- model.matrix(~0+group)
colnames(des) <- levels(group)
contr <- makeContrasts(a-b, levels = des) # comparison of control with drug response of susceptible strain of species1
# First set of 20 genes are genuinely differentially expressed
index1 <- 1:20
y[index1,4:6] <- y[index1,4:6]+1
# Second set of 20 genes are not DE
index2 <- 21:40
camera(y, list(set1=index1,set2=index2), design, contr = contr)
However, I am interested to know the interaction effect of species and drug. For example, if I am to compare gene expression profile of species1 control with species2 susceptible, then the difference in gene expression is partly because of drug treatment and partly because of being different species. How can I make the design table so as to get gene level scores (e.g. moderated t-statistics) that represent the drug*species interaction effect and then run camera?
I hope my problem is understandable!
Thanks in advance!
That was very helpful! Thanks a lot!!