Hello, I am working through the vignettes and examples for DESeq2 ?results
. I had a question(s) regarding lfcshrinkage. Looking at Example 2:
## Example 2: two conditions, two genotypes, with an interaction term
dds <- makeExampleDESeqDataSet(n=100,m=12)
dds$genotype <- factor(rep(rep(c("I","II"),each=3),2))
design(dds) <- ~ genotype + condition + genotype:condition
dds <- DESeq(dds)
resultsNames(dds)
# Note: design with interactions terms by default have betaPrior=FALSE
# the condition effect for genotype I (the main effect)
results(dds, contrast=c("condition","B","A"))
# the condition effect for genotype II
# this is, by definition, the main effect *plus* the interaction term
# (the extra condition effect in genotype II compared to genotype I).
results(dds, list( c("condition_B_vs_A","genotypeII.conditionB") ))
# the interaction term, answering: is the condition effect *different* across genotypes?
results(dds, name="genotypeII.conditionB")
How would you use the lfcShrink()
function in those three cases (in particular with 'apeglm', which won't work with contrast in the second results()
call)? Are you supposed to apply shrinkage to the interaction term? How would you specify a different alpha (0.05)?
Thanks! I saw the refactoring and I think I understand it, but I am curious about interpreting the outcome. I can for instance relevel
genotype
onII
and then it will become a singlecoef
term, but then would thecondition effect for this re-leveled genotype II
still be the sum of thecondition effect for genotype I
+interaction term
in the original model?For context, we are looking at a situation where a
treatment
had an expected behavioral effect only in males. So ultimately the goal is to look at DE genes in the males, but separate that gene set into groups for which the condition effect is different across sex, and those for which it is not. So with a~sex + treatment + sex:treatment
model, that would be"treatment_B_vs_A"
DE genes, separated by whether they are also DE in"sexF.treatmentB"
.Correct? or should I be refactoring this into a single
grp
variable and looking at something likegrp_Mtrt_vs_Muntrt
to get my DE gene set and then splitting them by whether they also are in the DE gene setgrp_Ftrt_vs_Mtrt
?This sound like the refactoring would be possible but introduce lots of placing where coding mistakes could happen. Why not just use ashr? It performs very well.