Gene set analysis with CAMERA with interaction effect
1
0
Entering edit mode
Suraj • 0
@suraj-7569
Last seen 9.7 years ago
Finland

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!

gene set analysis • 1.5k views
ADD COMMENT
1
Entering edit mode
@james-w-macdonald-5106
Last seen 12 hours ago
United States

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 COMMENT
0
Entering edit mode

That was very helpful! Thanks a lot!!
 

ADD REPLY

Login before adding your answer.

Traffic: 541 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6