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)