I have some questions about running a likelihood ratio test in the DESeq2 package in R. My experiment has two factors: sex and population, with 8-12 replicates per condition. I wish to identify genes that show an interaction between sex and population (specifically, genes with a sexually dimorphic expression in populations A and B, but not populations C and D).
Unfortunately, I'm not sure that I've set up the code properly, as my LRT results don't seem to corroborate with the results from some of my more heuristic methods.
Can anyone tell me if I'm making any glaring mistakes in my code?
countData <- read.table ("Counttable.txt", header=TRUE)
stickleDesign = data.frame (
row.names = colnames(countData),
sex = c("M", "F", "M", cont... )
pop = c("A", "A", "B", cont... )
dds <- DESeqDataSetFromMatrix (countData=countData,
design = ~sex)
dds <- DESeq(dds)
# Multi-Factor Analysis
design(dds) <- formula(~ sex + pop + sex:pop)
dds <- DESeq(dds, test=c("LRT"), full=~ sex + pop + sex:pop, reduced=~sex + pop)
resMFType <- results(dds, contrast=c("sex","M","F"))
Additionally, how do I tell that the model with the interaction term is necessarily better than one without? Is there an objective way to tell? I guess I was expecting a single p-value to represent the effectiveness of the model, but I don't see that I'm provided with anything like that.
Thanks for the help!