Hi all,
I have a set of RNA-seq data. It's a time series with 4 stages of leaf development, 2 treatments, and 3 bio-reps per stage x treat combo. I've already found DEGs for the treatments, and DEGs for the stage:treatment interaction, but what I would like to do is find DEGs at each stage comparison (e.g. genes that are specifically differentially expressed at stage 2, etc.). I've tried doing clustering with maSigPro, but it did not really provide the resolution that I need.
Any ideas on this? My code is below.
Thanks in advance!
reads <- as.matrix(read.csv("genereads_ONLY3.txt", sep = '\t', row.names = 1, header = TRUE))
meta <- read.table("metatest2.txt", header = TRUE)
mm <- model.matrix(~ stage + treat + stage:treat, data = meta)
mm2 <- model.matrix(~ stage + treat, data = meta)
mm3 <- model.matrix(~ stage, data = meta)
ddsd <- DESeqDataSetFromMatrix(countData = reads,
colData = meta,
design = mm)
dds <- ddsd[ rowSums(counts(ddsd)) > 10, ]
#running DESeq
dds.interaction <- DESeq(dds, test = "LRT", full = mm, reduced = mm2)
#Retrieves genes that are DEGs due to interaction of stage and treatment
dds.interact.results <- results(dds.interaction)
sum(dds.interact.results$padj < 0.05, na.rm = TRUE)
sub.interact <- subset(dds.interact.results, padj<0.05)
dds.interact.DEGs.genes <- row.names(sub.interact)
dds.interact.DEGs.genes.df <- as.data.frame(dds.interact.DEGs.genes)
colnames(dds.interact.DEGs.genes.df) <- "GENES"
dds.treatonly <- DESeq(dds, test = "LRT", full = mm, reduced = mm3)
dds.treat.DEGS <- results(dds.treatonly)
sum(dds.treat.DEGS$padj < 0.05, na.rm = TRUE)
sub.treat <- subset(dds.treat.DEGS, padj<0.05)
dds.treat.DEGs.genes <- row.names(sub.treat)
dds.treat.DEGs.genes.df <- as.data.frame(dds.treat.DEGs.genes)
colnames(dds.treat.DEGs.genes.df) <- "GENES"
DIFFERENCE <- setdiff(dds.treat.DEGs.genes.df, dds.interact.DEGs.genes.df)
Thank you for the answer! Whenever I use resultsNames(), I get the same thing for my different tests (the different tests have different designs, as can be seen above).
I get:
In which "X800" represents my treatment.
I know I've seen other outcomes of "results()" online that seem much lengthier. How would I go about designing my 'contrast = " with that output to get a comparison of a single stage between my control and treat?
Here is my full model matrix in case that helps:
When you run
DESeqDataSetFrom...()
it should be printing a message to the console that you have encoded the variables as numeric, and that you may have intended them to be factors instead. You should heed that message.(Edit: the message comes upon dataset construction, not
DESeq()
)Thank you! For some reason, that message never popped up. But that was in fact the issue!
Here is what happens when you create a
dds
with an integer-valued numeric in the design:What version of DESeq2 are you using?