I have a set of high-troughput experiments with 2 genotypes ("WT" and "prg1") and 3 treatments ("RT", "HS30" and "HS30RT120"), and there are 2 replicates for each of the genotype x treatment combinations.
The read counts for the genes are summarized in a file that I load as follows in R:
> counts_data <- read.table("path/to/my/file", header=TRUE, row.names="gene")
> colnames(counts_data)
[1] "WT_RT_1" "WT_HS30_1" "WT_HS30RT120_1" "prg1_RT_1"
[5] "prg1_HS30_1" "prg1_HS30RT120_1" "WT_RT_2" "WT_HS30_2"
[9] "WT_HS30RT120_2" "prg1_RT_2" "prg1_HS30_2" "prg1_HS30RT120_2"
I describe the experiments as follows:
> col_data <- DataFrame(
geno = c(rep("WT", times=3), rep("prg1", times=3), rep("WT", times=3), rep("prg1", times=3)),
treat = rep(c("RT", "HS30", "HS30RT120"), times=4),
rep = c(rep("1", times=6), rep("2", times=6)),
row.names = colnames(counts_data))
> col_data
DataFrame with 12 rows and 3 columns
geno treat rep
<character> <character> <character>
WT_RT_1 WT RT 1
WT_HS30_1 WT HS30 1
WT_HS30RT120_1 WT HS30RT120 1
prg1_RT_1 prg1 RT 1
prg1_HS30_1 prg1 HS30 1
... ... ... ...
WT_HS30_2 WT HS30 2
WT_HS30RT120_2 WT HS30RT120 2
prg1_RT_2 prg1 RT 2
prg1_HS30_2 prg1 HS30 2
prg1_HS30RT120_2 prg1 HS30RT120 2
I want to build a DESeq2 object that I could use to either:
- find differentially expressed genes when the treatment varies for a given fixed genotype
or:
- find differentially expressed genes when the genotype varies for a given fixed treatment
In the bioconductor help forum I think I've found a somewhat similar situation, and I read the following:
Try a design of ~ genotype + genotype:condition
Then you will have a condition effect for each level of genotype, including the reference level.
You can constrast pairs of them using the list style of the 'contrast' argument.
However, this doesn't explain how to apply this "list style" to the "contrast" argument. And the above situation seems to be asymmetrical. By that I mean that genotype and condition do not seem to have an interchangeable role.
So I tried the following more symmetric formula:
> dds <- DESeqDataSetFromMatrix(
countData = counts_data,
colData = col_data,
design = ~ geno + treat + geno:treat)
> dds <- DESeq(dds)
Now, can I for instance get the differential expression results when comparing treatment "HS30" against "RT" as a reference, in genotype "prg1"?
And how?
If I understand correctly, the above-mentioned "list style" uses names given by the resultsNames
function. In my case, I have the following:
> resultsNames(dds)
[1] "Intercept" "geno_WT_vs_prg1"
[3] "treat_HS30RT120_vs_HS30" "treat_RT_vs_HS30"
[5] "genoWT.treatHS30RT120" "genoWT.treatRT"
I guess I would need a contrast between "genoprg1.treatRT" and a "genoprg1.treatHS30", but these are not in the above results names.
I'm lost.