Question: Contrast for DESeq2 design with interaction and blocking factor
0
gravatar for erwan.scaon
4 months ago by
erwan.scaon10
erwan.scaon10 wrote:

Dear community,

I am trying to perform a differential abundance analysis on 16S data using DESeq2.

The metadata is as follow :
"Prelevement" : 2 types (P1, P2)
"Extraction" : 3 types (E1, E2, E3)
"NS" : Corresponds to donnor ID. Each donnor (10) had 6 samples, one for each case (P1 x E1, P1 x E2, P1 x E3, P2 x E1, P2 x E2, P2 x E3)

Thus design should be :
Extraction + Prelevement + Extraction:Prelevement + NS as blocking factor
design <- ~ NS + Extraction * Prelevement

I did read the ressources below with care : https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#interactions http://master.bioconductor.org/packages/release/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html ?results

But I am still unsure about how to extract specific contrasts from my dds. We want to retrieve differentially abundant taxa for :

# Prelevement effect :

P1 vs P2 for E1
P1 vs P2 for E2
P1 vs P2 for E3

# Extraction effect :

E1 vs E2 for P1
E1 vs E2 for P2
E1 vs E3 for P1
E1 vs E3 for P2
E2 vs E3 for P1
E2 vs E3 for P2

Ps : I did not set any reference level for my factors, given that there is no real "biologically relevant" reference here.
I am thus using R alphabetical order factors levels.

> levels(dds$Extraction)
[1] "E1"   "E2" "E3"
# "E1" as reference level
> levels(dds$Prelevement)
[1] "P1" "P2" 
# "P1" as reference level

Analysis was run with the following commands :

design <- ~ NS + Extraction * Prelevement
dds <- DESeqDataSetFromMatrix(countData = data, colData = meta, design = design)
sizeFactors(dds) <- res$GMPR
dds <- estimateDispersions(dds, fitType = "local")
dds <- nbinomWaldTest(dds, betaPrior = FALSE)

resultsNames(dds)

"Intercept"
" NSNS02vsNS01"
"NS
NS03vsNS01"
"NSNS04vsNS01"
"NS
NS05vsNS01"
"NSNS06vsNS01"
"NS
NS07vsNS01"
"NSNS08vsNS01"
"NS
NS09vsNS01"
"NSNS10vsNS01"
"Extraction
E2vsE1"
"ExtractionE3vsE1"
"Prelevement
P2vsP1"
"ExtractionE2.PrelevementP2"
"ExtractionE3.PrelevementP2"

Below are attempts to extract specific contrasts, feedback is welcomed :

## Prelevement effect :
# [P1 vs P2 for E1] for Extraction reference level (E1)
results(dds, contrast = c("Prelevement", "P1", "P2"))
# [P1 vs P2 for E2] for Extraction "E2" (main effect *plus* the interaction term)
results(dds, list(c("Prelevement_P1_vs_P2", "ExtractionE2.PrelevementP2")))
# [P1 vs P2 for E3] for Extraction "E3" (main effect *plus* the interaction term)
results(dds, list(c("Prelevement_P1_vs_P2", "ExtractionE3.PrelevementP2")))

## Extraction effect :
# [E1 vs E2 for P1]
results(dds, contrast = c("Extraction", "E1", "E2"))
# [E1 vs E3 for P1]
results(dds, contrast = c("Extraction", "E1", "E3"))
# [E2 vs E3 for P1]
results(dds, contrast = c("Extraction", "E2", "E3"))
# [E1 vs E2 for P2]
results(dds, name="ExtractionE2.PrelevementP2")
# [E1 vs E3 for P2]
?
# [E2 vs E3 for P2]
results(dds, contrast=list("ExtractionE2.PrelevementP2", "ExtractionE3.PrelevementP2"))

Sorry if some of these questions are trivial or if I am not always using the right terminology,
Best regards

deseq2 • 117 views
ADD COMMENTlink modified 4 months ago • written 4 months ago by erwan.scaon10
Answer: Contrast for DESeq2 design with interaction and blocking factor
1
gravatar for Michael Love
4 months ago by
Michael Love25k
United States
Michael Love25k wrote:

It looks like your desired contrasts could be made by combining all the factors, as described in the vignette?

ADD COMMENTlink written 4 months ago by Michael Love25k

Thanks for your reply,

Using DESeq2 version 1.22.2

meta$group <- factor(paste0(meta$Extraction, meta$Prelevement))
design <- ~ NS + group
dds <- DESeqDataSetFromMatrix(countData = data, colData = meta, design = design)
sizeFactors(dds) <- res$GMPR
dds <- estimateDispersions(dds, fitType = "local")
dds <- nbinomWaldTest(dds, betaPrior = TRUE)

dds$group

[1] E1P2 E1P2 E1P2 E1P2 E1P2 E1P2 E1P2 E1P2 E1P2 E1P2 E1P1 E1P1 E1P1 E1P1 E1P1 E1P1 E1P1 E1P1 E1P1 E1P1 E3P2 E3P2 E3P2 E3P2 E3P2 E3P2 E3P2 E3P2 E3P2 E3P2 E3P1
[32] E3P1 E3P1 E3P1 E3P1 E3P1 E3P1 E3P1 E3P1 E3P1 E2P2 E2P2 E2P2 E2P2 E2P2 E2P2 E2P2 E2P2 E2P2 E2P2 E2P1 E2P1 E2P1 E2P1 E2P1 E2P1 E2P1 E2P1 E2P1 E2P1
Levels: E1P1 E1P2 E2P1 E2P2 E3P1 E3P2

resultsNames(dds)

[1] "Intercept" "NSNS01" "NSNS02" "NSNS03" "NSNS04"
"NSNS05" "NSNS06" "NSNS07" "NSNS08" "NSNS09" "NSNS10"
"groupE1P1" "groupE1P2"
[14] "groupE2P1" "groupE2P2" "groupE3P1" "groupE3P2"

Notice that there is no more coefficients for interactions (with design ~ NS + group as compared to design ~ NS + Extraction * Prelevement)


Following your advice, I wrote these contrasts :

# Prelevement effect for Extraction E1
results(dds, contrast=c("group", "E1P2", "E1P1"))
# Prelevement effect for Extraction E2
results(dds, contrast=c("group", "E2P2", "E2P1"))
# Prelevement effect for Extraction E3
results(dds, contrast=c("group", "E3P2", "E3P1"))
# Extraction effect (E1 vs E2) for Prelevement P1
results(dds, contrast=c("group", "E1P1", "E2P1"))
# Extraction effect (E1 vs E2) for Prelevement P2
results(dds, contrast=c("group", "E1P2", "E2P2"))
# Extraction effect (E1 vs E3) for Prelevement P1
results(dds, contrast=c("group", "E1P1", "E3P1"))
# Extraction effect (E1 vs E3) for Prelevement P2
results(dds, contrast=c("group", "E1P2", "E3P2"))
# Extraction effect (E2 vs E3) for Prelevement P1
results(dds, contrast=c("group", "E2P1", "E3P1"))
# Extraction effect (E2 vs E3) for Prelevement P2
results(dds, contrast=c("group", "E2P2", "E3P2"))

Does this seems right ? Best regards

ADD REPLYlink modified 4 months ago • written 4 months ago by erwan.scaon10
1

Yes, aren’t those the group comparisons you wanted?

ADD REPLYlink written 4 months ago by Michael Love25k
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: 202 users visited in the last hour