Question

How to make design matrices or compare contrast in limma

0

Entering edit mode

mohammedtoufiq91 ▴ 10

@mohammedtoufiq91-17679

Last seen 5 weeks ago

United States

Hi, I am working with the transcriptomics expression data and trying to identify differentially expressed genes between "Molecules" groups thorugh a limma R package.

Currently, I have tried to perform simple comparison within the "Molecules" groups specific to "Condition" and "Stage" columns. For instance, I have extracted the samples of "Healthy" from "Condition" and "Stage" columns and compared between LPS vs. 6h (Baseline), and LTA_vs_6h (Baseline). Likewise for the "Group_A" and "Rem" sample-set. Both these designs are intra group comparisons. The sample metadata and example R code is as follows:

I am interested in inter-group comparisons; My question is, how can I make the contrasts or compare contrasts by mixing both the "Healthy" and "Group_A" sample-set simultaneously by considering both "6h samples" and "Healthy samples" as the Baseline.

Thank you, Toufiq

library(limma)
dput(Sample_metadata)
structure(list(Subject = c("HC1", "HC1", "HC1", "HC2", "HC2", 
                           "HC2", "P1", "P1", "P1", "P2", "P2", "P2", "P3", "P3", "P3"), 
               Molecules = c("6h", "LTA", "LPS", "6h", "LTA", "LPS", "6h", 
                             "LTA", "LPS", "6h", "LTA", "LPS", "6h", "LTA", "LPS"), Condition = c("Healthy", 
                                                                                                  "Healthy", "Healthy", "Healthy", "Healthy", "Healthy", "Group_A", 
                                                                                                  "Group_A", "Group_A", "Group_A", "Group_A", "Group_A", "Group_A", 
                                                                                                  "Group_A", "Group_A"), Stage = c("Healthy", "Healthy", "Healthy", 
                                                                                                                                   "Healthy", "Healthy", "Healthy", "Rem", "Rem", "Rem", "Rem", 
                                                                                                                                   "Rem", "Rem", "Rem", "Rem", "Rem")), class = "data.frame", row.names = c("S1", 
                                                                                                                                                                                                            "S2", "S3", "S4", "S5", "S6", "S7", "S8", "S9", "S10", "S11", 
                                                                                                                                                                                                            "S12", "S13", "S14", "S15"))
#>     Subject Molecules Condition   Stage
#> S1      HC1        6h   Healthy Healthy
#> S2      HC1       LTA   Healthy Healthy
#> S3      HC1       LPS   Healthy Healthy
#> S4      HC2        6h   Healthy Healthy
#> S5      HC2       LTA   Healthy Healthy
#> S6      HC2       LPS   Healthy Healthy
#> S7       P1        6h   Group_A     Rem
#> S8       P1       LTA   Group_A     Rem
#> S9       P1       LPS   Group_A     Rem
#> S10      P2        6h   Group_A     Rem
#> S11      P2       LTA   Group_A     Rem
#> S12      P2       LPS   Group_A     Rem
#> S13      P3        6h   Group_A     Rem
#> S14      P3       LTA   Group_A     Rem
#> S15      P3       LPS   Group_A     Rem

^{Created on 2022-04-10 by the [reprex package](https://reprex.tidyverse.org) (v2.0.1)}

Molecules <- factor(Sample_metadata$Molecules)

# simple design; considering samples from the "Healthy" condition only
design_Molecules_Healthy <- model.matrix(~ 0 + Molecules)
dput(design_Molecules_Healthy)
structure(c(1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 
            1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0), .Dim = c(12L, 
                                                                               3L), .Dimnames = list(c("1", "2", "3", "4", "5", "6", "7", "8", 
                                                                                                       "9", "10", "11", "12"), c("Molecules6h", "MoleculesLPS", "MoleculesLTA"
                                                                                                       )), assign = c(1L, 1L, 1L), contrasts = list(Molecules = "contr.treatment"))
#>    Molecules6h MoleculesLPS MoleculesLTA
#> 1            1            0            0
#> 2            0            0            1
#> 3            0            1            0
#> 4            1            0            0
#> 5            0            0            1
#> 6            0            1            0
#> 7            1            0            0
#> 8            0            0            1
#> 9            0            1            0
#> 10           1            0            0
#> 11           0            0            1
#> 12           0            1            0
#> attr(,"assign")
#> [1] 1 1 1
#> attr(,"contrasts")
#> attr(,"contrasts")$Molecules
#> [1] "contr.treatment"


## Make contrast
Cont.matrix_Molecules_Healthy <- makeContrasts("LPS_vs_0h" = MoleculesLPS-Molecules6h,  
                                                "LTA_vs_6h" = MoleculesLTA-Molecules6h,levels=design_Molecules_Healthy) 

dput(Cont.matrix_Molecules_Healthy)

structure(c(-1, 1, 0, -1, 0, 1), .Dim = 3:2, .Dimnames = list(
  Levels = c("Molecules6h", "MoleculesLPS", "MoleculesLTA"), 
  Contrasts = c("LPS_vs_0h", "LTA_vs_6h")))
#>               Contrasts
#> Levels         LPS_vs_0h LTA_vs_6h
#>   Molecules6h         -1        -1
#>   MoleculesLPS         1         0
#>   MoleculesLTA         0         1


# simple design; considering samples from the "Group A" condition + "Rem" stage only
design_Molecules_Group_A <- model.matrix(~ 0 + Molecules)
dput(design_Molecules_Group_A)

structure(c(1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 
            1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0), .Dim = c(12L, 
                                                                               3L), .Dimnames = list(c("1", "2", "3", "4", "5", "6", "7", "8", 
                                                                                                       "9", "10", "11", "12"), c("Molecules6h", "MoleculesLPS", "MoleculesLTA"
                                                                                                       )), assign = c(1L, 1L, 1L), contrasts = list(Molecules = "contr.treatment"))
#>    Molecules6h MoleculesLPS MoleculesLTA
#> 1            1            0            0
#> 2            0            0            1
#> 3            0            1            0
#> 4            1            0            0
#> 5            0            0            1
#> 6            0            1            0
#> 7            1            0            0
#> 8            0            0            1
#> 9            0            1            0
#> 10           1            0            0
#> 11           0            0            1
#> 12           0            1            0
#> attr(,"assign")
#> [1] 1 1 1
#> attr(,"contrasts")
#> attr(,"contrasts")$Molecules
#> [1] "contr.treatment"

## Make contrast
Cont.matrix_Molecules_Group_A <- makeContrasts("LPS_vs_0h" = MoleculesLPS-Molecules6h,  
                                               "LTA_vs_6h" = MoleculesLTA-Molecules6h,levels=design_Molecules_Group_A)
dput(Cont.matrix_Molecules_Group_A)
structure(c(-1, 1, 0, -1, 0, 1), .Dim = 3:2, .Dimnames = list(
  Levels = c("Molecules6h", "MoleculesLPS", "MoleculesLTA"), 
  Contrasts = c("LPS_vs_0h", "LTA_vs_6h")))
#>               Contrasts
#> Levels         LPS_vs_0h LTA_vs_6h
#>   Molecules6h         -1        -1
#>   MoleculesLPS         1         0
#>   MoleculesLTA         0         1

^{Created on 2022-04-10 by the [reprex package](https://reprex.tidyverse.org) (v2.0.1)}

designmatrices differentialexpression makeContrasts limma model.matrix • 1.3k views

ADD COMMENT • link 2.0 years ago mohammedtoufiq91 ▴ 10

1

Entering edit mode

Do you have multiple samples from the same subjects or is every sample from a different subject? How many independent subjects are involved in your study?

ADD REPLY • link 2.0 years ago Gordon Smyth 50k

0

Entering edit mode

Hi Gordon Smyth Thank you for the prompt response. I have 10 independent Healthy Controls and 25 independent Subjects. I have multiple samples (molecules) from the different subjects. I have revised the sample metadata above (please refer).

Thank you, Toufiq

ADD REPLY • link 2.0 years ago mohammedtoufiq91 ▴ 10

score 2 · Accepted Answer · 2022-04-10

You appear to have a multi-level experiment with multiple samples from each subject, so you should follow the section on multi-level experiments in the limma User's Guide.

You should not subset the data and you must account for correlation between repeat samples from the same subject.

I hope you don't mind if I give you some advice about the code presentation in your question. I see that you've tried to do the right thing by constructing a reproducible example using reprex, but in truth dput and reprex are not necessary in the context of this forum. It would sufficient and easier to us to read if you just print the data.frame to the terminal in the usual R way. Using dput for objects that are automatically created by the code is especially unnecessary. Also, we ideally need to see all the lines of the metadata instead of just a subset.