I'm using DESeq2 to do differentially expression analysis. I have 8 liver slices from control group, and 6 liver slices from treatment group. So the purpose is to find out the differentially expressed genes between these two groups. But those slices are not independent: slices from control group the slices from treatment group may come from the same fish.
I'm using DESeq2 package in R, so my colData looks like:
samples fishGroup expGroup 1 Sample12 3 control 2 Sample18 4 control 3 Sample25 5 control 4 Sample31 6 control 5 Sample44 8 control 6 Sample6 2 control 7 Sample11 2 treat 8 Sample17 3 treat 9 Sample24 4 treat 10 Sample30 5 treat 11 Sample37 6 treat 12 Sample50 8 treat
If I perform DESeq2 based on independent sample test, then only "expGroup" will be used:
dds_independent <- DESeqDataSetFromMatrix(countData, colData, design = ~expGroup)
But that will miss the samples' pair information. So I guess that the paired sample test is better.
In the section of "Model matrix not full rank" from
vignette('DESeq2'), it seems to discuss about this problem, to use a design like this:
design = ~ expGroup + fishGroup. I'm no expert on statistics, but it seems that the two factors are treated equally, but in my understanding, expGroup should be prior, because fishGroup only gives the pair information, since I want to find the differentially expressed genes in different experiment groups (control and treat).
I ran both (
design = ~expGroupand
design = ~expGroup + fishGroup) again, but But the result is quite abnormal: I got 10 differentially expressed (DE) genes (padj < 0.05) from the first design, but 2412 DE genes from the second design. I've looked into the DE genes from the second design, and they're not differentially expressed, since you'll find that the slice sample pairs (from the same fish) give you different conclusion about whether the gene is up regulated or down regulated. Take one gene which is regarded as DE gene from the second design:
fishGroup expGroup count 2 control 5.1 2 treat 25.6 3 control 222.9 3 treat 143.4 4 control 25.6 4 treat 164.0 5 control 0 5 treat 5.1 6 control 151.1 6 treat 61.5 8 control 886.3 8 treat 141
So anyone has any idea of how to design the "design" to tell DESeq2 to do a paired sample test?