Contrast for DESeq2 design with interaction and blocking factor
1
0
Entering edit mode
erwan.scaon ▴ 10
@erwanscaon-20285
Last seen 3.0 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"
"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 • 23k views
ADD COMMENT
1
Entering edit mode
@mikelove
Last seen 9 hours ago
United States

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

ADD COMMENT
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
1
Entering edit mode

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

ADD REPLY

Login before adding your answer.

Traffic: 969 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6