Question: Creating contrasts which do not include the reference level in limma, edgeR or DESeq2
0
gravatar for Ben
6 weeks ago by
Ben0
United States
Ben0 wrote:

Hi!

This question is regarding the creation of a contrast matrix in limma, edgeR or DESeq2, while analyzing RNA-Seq data. The question is: "How to create contrasts which represent comparisons between two factors which are not the reference level (intercept) using a model which contains an interaction term?"

(I do believe the interaction term is the important component here, since the interaction term will lead to a more complicated contrast matrix.)

Let's have a look at an example:

# Create data frame with information.
groups <- data.frame(group = c(rep(1,6), rep(2,6), rep(3,6), rep(4,6)),
                     subject = c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6,7,7,7,8,8,8),
                     treat = c(rep("a",6),rep("b",6),rep("a",6),rep("b",6)),
                     time = c(rep(c(1,2,3), 4), rep(c(1,4,5), 4)),
                     outcome = c(3,2,4,5,3,4,9,8,7,8,8,10,2,3,1,1,3,2,9,9,7,6,8,6))

# Setting levels for fixed effects.
groups$treat <- factor(groups$treat)
groups$time <- factor(groups$time)
groups$group <- factor(groups$group)
groups$subject <- factor(groups$subject)

# Creating design formula.
design <- model.matrix(outcome ~ treat * time, data = groups)
design
#>    (Intercept) treatb time2 time3 time4 time5 treatb:time2 treatb:time3
#> 1            1      0     0     0     0     0            0            0
#> 2            1      0     1     0     0     0            0            0
#> 3            1      0     0     1     0     0            0            0
#> 4            1      0     0     0     0     0            0            0
#> 5            1      0     1     0     0     0            0            0
#> 6            1      0     0     1     0     0            0            0
#> 7            1      1     0     0     0     0            0            0
#> 8            1      1     1     0     0     0            1            0
#> 9            1      1     0     1     0     0            0            1
#> 10           1      1     0     0     0     0            0            0
#> 11           1      1     1     0     0     0            1            0
#> 12           1      1     0     1     0     0            0            1
#> 13           1      0     0     0     0     0            0            0
#> 14           1      0     0     0     1     0            0            0
#> 15           1      0     0     0     0     1            0            0
#> 16           1      0     0     0     0     0            0            0
#> 17           1      0     0     0     1     0            0            0
#> 18           1      0     0     0     0     1            0            0
#> 19           1      1     0     0     0     0            0            0
#> 20           1      1     0     0     1     0            0            0
#> 21           1      1     0     0     0     1            0            0
#> 22           1      1     0     0     0     0            0            0
#> 23           1      1     0     0     1     0            0            0
#> 24           1      1     0     0     0     1            0            0
#>    treatb:time4 treatb:time5
#> 1             0            0
#> 2             0            0
#> 3             0            0
#> 4             0            0
#> 5             0            0
#> 6             0            0
#> 7             0            0
#> 8             0            0
#> 9             0            0
#> 10            0            0
#> 11            0            0
#> 12            0            0
#> 13            0            0
#> 14            0            0
#> 15            0            0
#> 16            0            0
#> 17            0            0
#> 18            0            0
#> 19            0            0
#> 20            1            0
#> 21            0            1
#> 22            0            0
#> 23            1            0
#> 24            0            1
#> attr(,"assign")
#>  [1] 0 1 2 2 2 2 3 3 3 3
#> attr(,"contrasts")
#> attr(,"contrasts")$treat
#> [1] "contr.treatment"
#> 
#> attr(,"contrasts")$time
#> [1] "contr.treatment"

Created on 2019-07-08 by the reprex package (v0.3.0)

This example shows a table containing the factor levels for the design matrix. After that, you can see the design formula and the design matrix printed to the console.

This dummy dataset contains an outcome (Y), which would be not the case for an RNA-Seq experiment, but I used this example for testing where I wanted an outcome (Y) - so bear with me here.

The intercept represents our reference level, which is "treat a" and "time 1". "treat b" (column 2 in the design matrix) represents the difference of "treat b / time 1" compared to "treat a / time 1". "time 2, 3, 4, 5" represent the respective time differences of "treat a".

Columns 7 to 10 in the design matrix represent the additional effects in "treat b".

Goal is to create a contrast which represents the comparison of "time 5" compared to "time 1" of "treat b". Here, "treat b / time 1" should be our reference level.

I do believe this should be the effect of "treat a / time 5" plus the additional effect captured within the coefficient "treatb:time5", so my contrast would look like this:

makeContrast(contrast = time5 + treatb:time5, ...)

This results in a contrast matrix which does not equals zero, which for me looks wrong.

I did find explanations about expressing "differences of differences" using a contrast matrix (e.g. here), but I do not want to express differences of differences. I do want the logFC and p-values for the comparison of "treat b / time 5" compared to "treat b / time 1".

Thanks for your help!

limma edger deseq2 • 108 views
ADD COMMENTlink modified 6 weeks ago by James W. MacDonald50k • written 6 weeks ago by Ben0

I think you need a design formula without intercept (~0 + ...) see edgeR manual (https://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf) p.69 section 4.4.6

ADD REPLYlink written 6 weeks ago by ATpoint10

Thanks for your answer and the reference!

Unfortunately, the example given in the edgeR document deals with model containing a single factor, not an interaction. Things are getting more complicated (at least from my perspective) with an interaction term. I will look into the option without an intercept, but I will keep this question open for now since I have the feeling the question isn't really answered in the context of interaction terms.

ADD REPLYlink written 6 weeks ago by Ben0
Answer: Creating contrasts which do not include the reference level in limma, edgeR or D
1
gravatar for James W. MacDonald
6 weeks ago by
United States
James W. MacDonald50k wrote:

So you are just asking the same question you already asked here, to which I said to combine the treatment and time. Evidently that didn't make any sense since you have just asked the same question again using a different post name (hint: don't do that - it's considered rude to just keep posting the same question in different posts).

Since words evidently didn't make sense, let's try an example.

> trttime <- paste(groups$treat, groups$time, sep = "_")
> design <- model.matrix(~0+trttime)
> colnames(design) <- gsub("trttime", "", colnames(design))
> design
   a_1 a_2 a_3 a_4 a_5 b_1 b_2 b_3 b_4 b_5
1    1   0   0   0   0   0   0   0   0   0
2    0   1   0   0   0   0   0   0   0   0
3    0   0   1   0   0   0   0   0   0   0
4    1   0   0   0   0   0   0   0   0   0
5    0   1   0   0   0   0   0   0   0   0
6    0   0   1   0   0   0   0   0   0   0
7    0   0   0   0   0   1   0   0   0   0
8    0   0   0   0   0   0   1   0   0   0
9    0   0   0   0   0   0   0   1   0   0
10   0   0   0   0   0   1   0   0   0   0
11   0   0   0   0   0   0   1   0   0   0
12   0   0   0   0   0   0   0   1   0   0
13   1   0   0   0   0   0   0   0   0   0
14   0   0   0   1   0   0   0   0   0   0
15   0   0   0   0   1   0   0   0   0   0
16   1   0   0   0   0   0   0   0   0   0
17   0   0   0   1   0   0   0   0   0   0
18   0   0   0   0   1   0   0   0   0   0
19   0   0   0   0   0   1   0   0   0   0
20   0   0   0   0   0   0   0   0   1   0
21   0   0   0   0   0   0   0   0   0   1
22   0   0   0   0   0   1   0   0   0   0
23   0   0   0   0   0   0   0   0   1   0
24   0   0   0   0   0   0   0   0   0   1
attr(,"assign")
 [1] 1 1 1 1 1 1 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$trttime
[1] "contr.treatment"

Now you can compare any group to any other group, without having to interpret what the different interaction terms mean. This obviously assumes that you don't care that A) you apparently have repeated measures from different subjects, and B) there is this outcome that you included but ignored as well.

ADD COMMENTlink written 6 weeks ago by James W. MacDonald50k

Thanks for your answer and help. For me - though - these two questions were not similar. Maybe this wasn't clear enough from my side, so I apologize. The question you are linking was regarding missing coefficients in the design model when using a model with no intercept.

The question here was about contrasts and setting contrasts in a model with intercept.

So, earlier question related to missing coefficients in model without intercept, this question regarding contrasts in model with intercept.

But again, maybe the way I formulated the question was not clear enough to show this difference.

Either way, you are making some correct assumptions, which are that I will do have a more complex model at the end. For demonstration purposes I chose a much simpler model, and yes, I also included an outcome here. I will take care of the repeated measures. That being said, simply removing my intercept in the answer you provided and merging my two fixed effects into one fixed effect was not the answer I was looking for, even though it is - in this demonstration example - a valid solution.

In the meantime, I found my answer regarding creating the contrast I want while maintaining the model as it is in my initial question.

Thanks much again!

ADD REPLYlink written 5 weeks ago by Ben0
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 16.09
Traffic: 218 users visited in the last hour