Help with model.matrix and creating the right contrast matrix for limma
1
0
Entering edit mode
rober204 • 0
@b48e8083
Last seen 12 days ago
Canada

Hi Limma experts, This is my first time using Limma and I am struggling to determine which contrasts to construct to answer my research questions/ if I have chosen the correct design matrix to most effectively answer my questions. The two questions I have:

  1. Is there a treatment effect?
  2. Is there a temporal effect?

My data and study design: Treatment (control, low, high), and Time (day minus 5, day 5, day 15, and day 30). The same pool of subjects are repeatedly sampled at each time point within each treatment, but they are sampled independently of each other (there is no way to know if it is the same individuals are being sampled at each time point). It's important to note that the experiment started on day 0, so day minus 5 (m5) is a temporal reference, prior to exposure.

I have decided to use model.matrix(~0 + group) after referencing Law (2020), and constructed groups: C_m5, C_5, C_15, C_30, L_m5, L_5, L_15, L_30, H_m5, H_5, H_15, H_30. However, there are a few other options I believe could be more useful for my case including model.matrix(~treatment*time) or model.matrix(~0+treatment+treatment:time), but I do not understand the workings of contrasts well enough to know which is best.

For my current model I chose (model.matrix(~0 + group)), I created a contrast matrix with all pairwise comparisons between groups... which looks pretty disturbing as you'll see below. I'm very unsure if this is the most effective way to make the contrasts, or if i would be better off making contrasts in another fashion, such as using interactions between treatment and time with one of the other models.

I'll Include the code of what I have thus far for limma, it's pretty ugly with all of the contrasts, which makes me feel it is wrong. I am hoping that someone is able to provide some feedback on what i have done, and if there is a different direction I should be taking. Further, with the results from decideTests, I am unable to make a vennDiagram, I will show the error below, but i am assuming it is because of how I made the contrasts, so it is not picking up my two factors (treatment and time).

Thank you in advanced for taking the time to help me today!

group <- interaction(treatment, time)

mm <- model.matrix(~0 + group) 
fit <- lmFit(log_mucus_cL, mm) 

#all pairwise comparisons...
cm <- makeContrasts(groupC._5vsgroupH._5 = groupH._5 - groupC._5,
                    groupC._5vsgroupL._5 = groupL._5 - groupC._5,
                    groupC._5vsgroupC.15 = groupC.15 - groupC._5,
                    groupC._5vsgroupH.15 = groupH.15 - groupC._5,
                    groupC._5vsgroupL.15 = groupL.15 - groupC._5,
                    groupC._5vsgroupC.30 = groupC.30 - groupC._5,
                    groupC._5vsgroupH.30 = groupH.30 - groupC._5,
                    groupC._5vsgroupL.30 = groupL.30 - groupC._5,
                    groupC._5vsgroupC.m5 = groupC._5 - groupC.m5,
                    groupC._5vsgroupH.m5 = groupC._5 - groupH.m5,
                    groupC._5vsgroupL.m5 = groupC._5 - groupL.m5,
                    groupH._5vsgroupL._5 = groupH._5 - groupL._5,
                    groupH._5vsgroupC.15 = groupH._5 - groupC.15,
                    groupH._5vsgroupH.15 = groupH.15 - groupH._5,
                    groupH._5vsgroupL.15 = groupH._5 - groupL.15,
                    groupH._5vsgroupC.30 = groupH._5 - groupC.30,
                    groupH._5vsgroupH.30 = groupH.30 - groupH._5,
                    groupH._5vsgroupL.30  = groupH._5 - groupL.30,
                    groupH._5vsgroupC.m5  = groupH._5 - groupC.m5,
                    groupH._5vsgroupH.m5  = groupH._5 - groupH.m5,
                    groupH._5vsgroupL.m5  = groupH._5 - groupL.m5,
                    groupL._5vsgroupC.15  = groupL._5 - groupC.15,
                    groupL._5vsgroupH.15  = groupH.15 - groupL._5,
                    groupL._5vsgroupL.15  = groupL.15 - groupL._5,
                    groupL._5vsgroupC.30  = groupL._5 - groupC.30,
                    groupL._5vsgroupH.30  = groupH.30 - groupL._5,
                    groupL._5vsgroupL.30  = groupL.30 - groupL._5,
                    groupL._5vsgroupC.m5  = groupL._5 - groupC.m5,
                    groupL._5vsgroupH.m5  = groupL._5 - groupH.m5,
                    groupL._5vsgroupL.m5  = groupL._5 - groupL.m5,
                    groupC.15vsgroupH.15  = groupC.15 - groupH.15,
                    groupC.15vsgroupL.15  = groupC.15 - groupL.15,
                    groupC.15vsgroupC.30  = groupC.30 - groupC.15,
                    groupC.15vsgroupH.30  = groupH.30 - groupC.15,
                    groupC.15vsgroupL.30  = groupL.30 - groupC.15,
                    groupC.15vsgroupC.m5  = groupC.15 - groupC.m5,
                    groupC.15vsgroupH.m5  = groupC.15 - groupH.m5,
                    groupC.15vsgroupL.m5  = groupC.15 - groupL.m5,
                    groupH.15vsgroupL.15  = groupH.15 - groupL.15,
                    groupH.15vsgroupC.30  = groupH.15 - groupC.30,
                    groupH.15vsgroupH.30  = groupH.30 - groupH.15,
                    groupH.15vsgroupL.30  = groupH.15 - groupL.30,
                    groupH.15vsgroupC.m5  = groupH.15 - groupC.m5,
                    groupH.15vsgroupH.m5  = groupH.15 - groupH.m5,
                    groupH.15vsgroupL.m5  = groupH.15 - groupL.m5,
                    groupL.15vsgroupC.30 =  groupL.15 - groupC.30,
                    groupL.15vsgroupH.30 =  groupH.30 - groupL.15,
                    groupL.15vsgroupL.30 =  groupL.30 - groupL.15,
                    groupL.15vsgroupC.m5 =  groupL.15 - groupC.m5,
                    groupL.15vsgroupH.m5 =  groupL.15 - groupH.m5,
                    groupL.15vsgroupL.m5 =  groupL.15 - groupL.m5,
                    groupC.30vsgroupH.30 =  groupH.30 - groupC.30,
                    groupC.30vsgroupL.30 =  groupL.30 - groupC.30,
                    groupC.30vsgroupC.m5 =  groupC.30 - groupC.m5,
                    groupC.30vsgroupH.m5 =  groupC.30 - groupH.m5,
                    groupC.30vsgroupL.m5 =  groupC.30 - groupL.m5,
                    groupH.30vsgroupL.30 =  groupH.30 - groupL.30,
                    groupH.30vsgroupC.m5 =  groupH.30 - groupC.m5,
                    groupH.30vsgroupH.m5 =  groupH.30 - groupH.m5,
                    groupH.30vsgroupL.m5 =  groupH.30 - groupL.m5,
                    groupL.30vsgroupC.m5 =  groupL.30 - groupC.m5,
                    groupL.30vsgroupH.m5 =  groupL.30 - groupH.m5,
                    groupL.30vsgroupL.m5 =  groupL.30 - groupL.m5,
                    groupC.m5vsgroupH.m5 =  groupC.m5 - groupH.m5,
                    groupC.m5vsgroupL.m5 =  groupC.m5 - groupL.m5,
                    groupH.m5vsgroupL.m5 =  groupH.m5 - groupL.m5,
                    levels = mm)

# Fit the contrasts
fit2 <- contrasts.fit(fit, contrasts = cm)

# Calculate the t-statistics for the contrasts
fit2 <- eBayes(fit2)

# Summarize results
results <- decideTests(fit2)
summary(results)

 library(VennDiagram)
vennDiagram(results)

error messages
#Error in matrix(0, noutcomes, ncontrasts) :
#invalid 'nrow' value (too large or NA)

sessionInfo( )
limma • 258 views
ADD COMMENT
0
Entering edit mode
@gordon-smyth
Last seen 2 hours ago
WEHI, Melbourne, Australia

Your question is not about the use of limma software or even about statistics. Your question is about how to conduct a scientific analysis of your specific experiment.

The purpose of forming a contrast and conducting a statistical test is to answer a scientific question. The problem with the contrasts that you have currently formed is that they have no relationship to the scientific questions that you have posed. Hence you will just get a mess of uninterpretable and possibly conflicting information. The reason why you got an error from vennDiagram is that there is no such thing as a Venn Diagram for 66 contrasts -- such a thing is just geometrically impossible.

If I was analysing an experiment like this, I would test for genes that vary over time for each treatment. That would lead to three tests, one for each treatment. Then I test for genes whose pattern over time varies between the treatments. That is one more test. That approach would involve four tests in total rather than the 66 tests you have now.

Or you could turn it around and instead test for treatment differences at each time. That would be four tests. It's up to you, depending on your scientific hypotheses and on the conclusions you want to draw from the data.

ADD COMMENT
0
Entering edit mode

Hi Gordon, thanks for taking time to respond to my question. Yep, I'm trying to do the latter for my analysis. I know how I want my scientific analysis to look, so my question is how to construct the contrasts to answer those questions, and if I was using the right model to do it effectively. As I said, I just started using this package, the knowledge gap here is in generating the contrasts to answer my question. In my mind, testing for treatment differences at each time would look like this, which is 12 tests instead of the 4 you stated, so clearly I am missing something about contrasts:

cm <- makeContrasts(
  #day -5
  HvsC_m5 = H.m5 - C.m5,
  HvsL_m5 = H.m5 - L.m5,
  LvsC_m5 = L.m5 - C.m5,

  #day 5
  HvsC_5 = H._5 - C._5,
  HvsL_5 = H._5 - L._5,
  LvsC_5 = L._5 - C._5,

  #day 15
  HvsC_15 = H.15 - C.15,
  HvsL_15 = H.15 - L.15,
  LvsC_15 = L.15 - C.15,

  #day 30 
  HvsC_30 = H.30 - C.30,
  HvsL_30 = H.30 - L.30,
  LvsC_30 = L.30 - C.30,
)

Further, this doesn't test for temporal differences within treatments, which would be further tests. Perhaps there is an example of writing contrasts in a similar matter that you can direct me to?

Again, thank you for your time.

ADD REPLY
0
Entering edit mode

To test for DE at day -5 you could use:

topTable(fit2, coef=c("LvsC_m5","HvsC_m5"))

That yields one anova test for each gene. You could choose any pair of contrasts at the same time to get the same test results.

There are many equivalent ways to setup the design matrix and contrasts for the same purpose, but this way is an easy as any and has the advantage of being explicit.

ADD REPLY

Login before adding your answer.

Traffic: 229 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