DESeq2 - desing, multiple groups - treatment in different time

Question

[LatifaP]

Hello, If someone can help I will be very grateful. My question is about which design should have my data to run deseq2 and see some comparisons.

First, I have the following data: 2genotypes, 2 tissue, and a treatment evaluated in 4 different times (control, t1, t2 and t3). In a matrix it is like this:

• library genotype tissue treatment replicate
• 1 A tissue1 Control 1
• 2 A tissue1 Time1 1
• 3 A tissue1 Time2 1
• 4 A tissue1 Time3 1
• 5 A tissue2 Control 1
• 6 A tissue2 Time1 1
• 7 A tissue2 Time2 1
• 8 A tissue2 Time3 1
• 9 A tissue1 Control 2
• 10 A tissue1 Time1 2
• 11 A tissue1 Time2 2
• 12 A tissue1 Time3 2
• 13 A tissue2 Control 2
• 14 A tissue2 Time1 2
• 15 A tissue2 Time2 2
• 16 A tissue2 Time3 2
• 17 A tissue1 Control 3
• 18 A tissue1 Time1 3
• 19 A tissue1 Time2 3
• 20 A tissue1 Time3 3
• 21 A tissue2 Control 3
• 22 A tissue2 Time1 3
• 23 A tissue2 Time2 3
• 24 A tissue2 Time3 3

I have the same for genotype B. So in total I have 48 samples.

I would like to know how gene expression changes over control, time1, time2, and time3 in one genotype and specific tissue, for example in the genotype B, in the tissue1.

So, I was applying this design:

dds <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = directory,
                                  design = ~genotype + tissue + treatment)
dds$group <- factor(paste0(dds$genotype,dds$tissue,dds$treatment))
design(dds) <- ~ group
dds <- DESeq(dds)  

then: resultsNames(dds):
 [1] "Intercept"                                                              "group_ATissue1T1_vs_ATissue1Control"     
 [3]  "group_ATissue1T2_vs_ATissue1Control"           "group_ATissue1T3_vs_ATissue1Control"     
 [5]  "group_ATissue2Control_vs_ATissue1Control"    "group_ATissue2T1_vs_ATissue1Control"     
 [7]  "group_ATissue2T2_vs_ATissue1Control"            "group_ATissue2T3_vs_ATissue1Control"     
 [9]  "group_BTissue1Control_vs_ATissue1Control"    "group_BTissue1T1_vs_ATissue1Control"     
[11] "group_BTissue1T2_vs_ATissue1Control"            "group_BTissue1T3_vs_ATissue1Control"     
[13] "group_BTissue2Control_vs_ATissue1Control"    "group_BTissue2T1_vs_ATissue1Control"     
[15] "group_BTissue2T2_vs_ATissue1Control"            "group_BTissue2T3_vs_ATissue1Control"

There I don´t have "group_BTissue1T1_vs_BTissue1Control", to compare the fold change between the T1 and the control in the tissue1 in the genotype B.

However, I do not know why, but I can run this:

res <- results(dds, alpha = 0.05, contrast = c("group", "BTissue1T1","BTissue1Control"))

But to generate MA plot, I need to do the shrinkage (lfcShrink), and for that, I need the coef :
coef = c("group_BTissue1T1_vs_BTissue1Control") which I don't have.

Please how should be my design to see the fold change across the different treatment, in a specific genotype and tissue?. I was employing the default wald test, but because treatment is evaluated in 4 different times (control, t1, t2 and t3), maybe I should have to employ the LRT test.

Thanks in advance for your time.

[LatifaP]

thanks!, but after doing the relevel then I can´t see other combinations

[ATpoint]

You need a design without an intercept to make all possible contrasts. I am not familiar enough with DESeq2 for hands-on help, but in general you need design = ~ 0 + (...). I am sure there is online help for this. In edgeR it is easy with the makeContrasts function, pretty sure there is something similar in DESeq2.

...and please comment to answers with ADD COMMENT, not the answer field, to keep things logically organized.

Answer 1

[swbarnes2]

See how all of your ResultsNames end in ATissue1Control?

Put this

dds$group <- relevel(dds$group, "BTissue1Control")

before

dds <- DESeq(dds)

[LatifaP]

Yes, thanks!, but after doing relevel, I can´t see other combinations.

Answer 2

[LatifaP]

I just decided to do another desing:

design = ~genotype * tissue * treatment

dss = DESeq(dss)

At the end of running I got this message: "10 rows did not converge in beta, labelled in mcols(object)$betaConv. Use larger maxit argument with nbinomWaldTest"

For resultsNames(dds), I got: [1] "Intercept" "genotype_B_vs_A" "tissue_Tissue2_vs_Tissue1"
[4] "treatment_T1_vs_Control" "treatment_T2_vs_Control" "treatment_T3_vs_Control"
[7] "genotypeB.tissueTissue2" "genotypeB.treatmentT1" "genotypeB.treatmentT2"
[10] "genotypeB.treatmentT3" "tissueTissue2.treatmentT1" "tissueTissue2.treatmentT2"
[13] "tissueTissue2.treatmentT3" "genotypeB.tissueTissue2.treatmentT1" "genotypeB.tissueTissue2.treatmentT2" [16] "genotypeB.tissueTissue2.treatmentT3"

I was trying to understand how to do the comparison, but I do not sure if I am understanding it well.

So, in the case of having: Base levels: A, tissue1 and Control ,

To observe:

T1 vs Control in the genotype A, specifically in tissue1

res <- results(dds, alpha = 0.05, contrast = c("treatment","T1","Control")) .............is it correrct?

T1 vs Control in the genotype A, specifically in tissue2

res <- results(dds, alpha = 0.05, contrast = list(c("treatment_T1_vs_Control","tissueTissue2.treatmentT1" ))) .............is it correrct?

T1 vs control in the genotype B, specifically in tissue1

res <- results(dds, alpha = 0.05, contrast = list(c("treatment_T1_vs_Control","genotypeB.treatmentT1" ))) .............is it correrct?

T1 vs control in the genotype B, specifically in tissue2

res <- results(dds, alpha = 0.05, contrast = list(c("treatment_T1_vs_Control","genotypeB.tissueTissue2.treatmentT1"))) .............is it correrct?

So, I have some problems:

• First, I dont sure If what I have done is well or not.

• Second, if I continue doing that I have to do several comparisons.

For this reason, I should employ the LRT test, but I do not know how my design should be and what should I put in reduce to see if there is any difference in the expression among the 4 different times of the treatment for an specific genotype and in an specific tissue.