How to interpret DESeq2 result?

Question

[helen]

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?

Answer 1

[Asaf]

1. You mean Sex, not Gender
2. Why are you expecting a lot of genes? Are M:B and F:B so different? More than M:A and F:A? Try to plot a PCA of the data and see the structure of it.
3. Getting only 13 genes is usually great news and it means you pinpointed the exact effect, getting a lot of genes means your conditions are too different to learn something useful.

[helen]

1. Yes, Sex would be more appropriate. Thanks for pointing it out. For 2. & 3. I think I misunderstood what the result is trying to tell. I thought this design is to compare condition A vs. B under specific sex, but in fact this interpretation is not correct. If I understand correctly, the fewer DEGs indicates the sex difference effect is very minimal in this experiment. 13 genes now seems an acceptable number.

[Asaf]

Yes, Sex:Condition interaction will give you the differences between A and B which are sex specific, in your case there are 13 genes that are differentially expressed between A and B in males only.