Question: How to interpret DESeq2 result?
1
5 weeks ago by
helen30
helen30 wrote:

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
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?

rna-seq deseq2 • 214 views
modified 16 days ago • written 5 weeks ago by helen30
3
5 weeks ago by
Asaf8.3k
Israel
Asaf8.3k wrote:
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.
1
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.