Search
Question: Gene set analysis with CAMERA with interaction effect
0
gravatar for Suraj
2.6 years ago by
Suraj0
Finland
Suraj0 wrote:

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!

ADD COMMENTlink modified 2.6 years ago by James W. MacDonald45k • written 2.6 years ago by Suraj0
1
gravatar for James W. MacDonald
2.6 years ago by
United States
James W. MacDonald45k wrote:

That's not really how you do GSEA. While that is the example from the help page for camera(), the example is an artificial construct that is simply intended to show what the results look like if your gene set is up-regulated or not.

The idea is to see if a set of genes (usually from some external source, like a KeGG pathway, or a GO term or whatever) is in aggregate being affected in your experiment. In other words, you would come up with a set of genes that you think are meaningful in the context of your experiment, and you are then testing to see if that set is in aggregate being affected by the treatment. One of the original rationales for this method was to allow comparisons between two different experiments, but that was in the context of two different labs making the same or similar comparisons, but having little overlap between their two sets of 'top genes'. If you have a set of arrays that you processed at the same time, and are normalizing together, then you don't need GSEA.

One thing you can do is just make the comparisons you care about, and then do Venn diagrams showing the genes that are consistent between two comparisons, and those that are unique. So for your question of interest you could do

contrast <- makeContrasts(a-b, d-e, levels = design)
fit <- lmFit(eset, design)
fit2 <- contrasts.fit(fit, contrast)
fit2 <- eBayes(fit2)
rslt <- decideTests(fit2)
vennDiagram(rslt)

The genes in the intersection will be those genes that are significant in both contrasts, and the genes outside the intersection are unique to each contrast. But note that the genes in the union are not required to have the same directionality (e.g., a gene that is up in a-b and a gene that is down in d-e will end up in the intersection).

You can also test the interaction, but that is a different question. The interaction tests for a difference in reaction to the drug that is dependent on species. You can certainly test that:

contrast <- makeContrasts(a-b-d+e, levels = design)
fit <- lmFit(eset, design)
fit2 <- contrasts.fit(fit, contrast)
fit2 <- eBayes(fit2)
topTable(fit2)

And you can always use camera() to see if a given gene set appears to be affected in aggregate, but in the context of an interaction this becomes much more difficult to interpret.

ADD COMMENTlink written 2.6 years ago by James W. MacDonald45k

That was very helpful! Thanks a lot!!
 

ADD REPLYlink written 2.6 years ago by Suraj0
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 314 users visited in the last hour