Hi all,

I have 6 subjects 1-6, and would like to identify the differentially expressed genes in condition A vs.B and consider the gender effect.

I arrange my coldata as follow to mimic the example in this section of DESeq2 manual: http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#group-specific-condition-effects-individuals-nested-within-groups , which gender can be approximate to grp, subject to ind, condition to cnd

gender  subject  condition
     M        1          A
     M        1          B
     M        2          A
     M        2          B
     M        3          A
     M        3          B
     F        4          A
     F        4          B
     F        5          A
     F        5          B
     F        6          A
     F        6          B

Add a column subject.n

coldata$subject.n <- factor(rep(rep(1:3,each=2),2))

so the new coldata would be:

gender  subject  condition  subject.n
     M        1          A          1
     M        1          B          1
     M        2          A          2
     M        2          B          2
     M        3          A          3
     M        3          B          3
     F        4          A          1
     F        4          B          1
     F        5          A          2
     F        5          B          2
     F        6          A          3
     F        6          B          3

and the columns of my model.matrix(~ gender + gender:subject.n + gender:condition, coldata) are

[1] "(Intercept)"             "gendermale"                  
[3] "genderfemale:subject.n2" "gendermale:subject.n2"    
[5] "genderfemale:subject.n3" "gendermale:subject.n3"    
[7] "genderfemale:conditionB" "gendermale:conditionB"

so my design formula would be

dds <- DESeqDataSetFromMatrix(countData = cts,
                              colData = coldata,
                              design = ~ gender + gender:subject.n + gender:condition)

then I filter out low read counts and perform differential expression analysis

keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
dds <- DESeq(dds)

and then extract the result

res <- results(dds, contrast=list("gendermale.conditionB","genderfemale.conditionB"))

However,

out of 18461 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 13, 0.07%
LFC < 0 (down)     : 0, 0%
outliers [1]       : 0, 0%
low counts [2]     : 11095, 60%
(mean count < 565)

only 13 up-regulated genes are found (I was expecting hundreds to thousands of DE genes)

Does anyone know how to interpret this result?

Does this mean under the gender-specific condition effect, comparing conditions A to B would render 13 up-regulated genes? Does this method reach my aim? And why there are so few DE genes?