Question

Contrast for DESeq2 design with interaction and blocking factor

0

Entering edit mode

erwan.scaon ▴ 10

@erwanscaon-20285

Last seen 3.7 years ago

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"
"NSNS03vsNS01"
"NSNS04vsNS01"
"NSNS05vsNS01"
"NSNS06vsNS01"
"NSNS07vsNS01"
"NSNS08vsNS01"
"NSNS09vsNS01"
"NSNS10vsNS01"
"ExtractionE2vsE1"
"ExtractionE3vsE1"
"PrelevementP2vsP1"
"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 • 25k views

ADD COMMENT • link 5.5 years ago erwan.scaon ▴ 10

score 1 · Answer 1 · 2019-06-11

1

Entering edit mode

Michael Love 43k

@mikelove

Last seen 1 day ago

United States

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

ADD COMMENT • link 5.5 years ago Michael Love 43k

0

Entering edit mode

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 REPLY • link 5.5 years ago erwan.scaon ▴ 10

1

Entering edit mode

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

ADD REPLY • link 5.5 years ago Michael Love 43k