Hi,

I have a quick question regarding plots. I have previously created a plot such as this to look at the log2fold change between treatments, but I now have an interaction effect between treatment and different time points. I attempted to recreate the same plot, 9 times, as to compare treatment T1 and TreatmentT2 etc and see how the plot changed but the plot is always the same. Therefore I am 100% sure I did this wrong.

How can I go about plotting the differences in my treatments over time? Would I need to do a heatmap? Can I actually do the graph as stated below in the code if so, how do I get it to properly display what I need it to.

Here is the code I have used and the graphing code for one of the timepoints. Please let me know how I can go about addressing this an or if there is another resource I can use to look at this.

dds = phyloseq_to_deseq2(physeq, design =  ~ Treatment + TimePoint + Treatment:TimePoint)

*#Set reference group for comparison------------------------------------------------------*
dds$Treatment<-relevel(dds$Treatment,ref="Unburned")
levels(dds$Treatment)

dds$TimePoint<-relevel(dds$TimePoint,ref="T1")
levels(dds$TimePoint)

***# calculate geometric means prior to estimate size factors-------required because data has zero's-----------------***
gm_mean = function(x, na.rm=TRUE){
  exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
}

geoMeans<-apply(counts(dds), 1, gm_mean)
dds<-estimateSizeFactors(dds, geoMeans = geoMeans)

***#Run DESeq and calculate values required-------------------------------------***
dds<-DESeq(dds, test= "LRT", fitType="local", reduced = ~ Treatment)
resultsNames(dds)

res = results(dds, name="Treatment-Burned_T2") 
res1 = res[order(res$padj, na.last=NA), ]
alpha = 0.05
sigtab = res[(res$padj < alpha), ]
sigtab = cbind(as(sigtab, "data.frame"), as(tax_table(physeq)[rownames(sigtab), ], "matrix"))
sigtab


****#This is the graph I am trying to create** 
**but when I look at the different contrast names I see the exact same results. Did I mess up above in the Deseq calculations or am I doing the contrast incorrectly?**

scale_fill_discrete <- function(palname = "Set1", ...) {
  scale_fill_brewer(palette = palette, ...)
}%>%
sigtabgen = subset(sigtab, !is.na(Phylum))+theme_set(theme_bw())

***# Phylum***
x<-tapply(sigtabgen$log2FoldChange, sigtabgen$Phylum, function(x) max(x))
x<-sort(x, TRUE)
sigtabgen$Phylum <- factor(as.character(sigtabgen$Phylum), levels = names(x))
***# Genus***
x<-tapply(sigtabgen$log2FoldChange, sigtabgen$Genus, function(x) max(x))
x<-sort(x, TRUE)
sigtabgen$Genus<-factor(as.character(sigtabgen$Genus), levels = names(x))

***#Plot----------------------------------------------------------------------------------------------------------------------------***
ggplot(sigtabgen, aes(x = Genus, y = log2FoldChange, color = Phylum)) + geom_point(size=6) + 
      geom_point(size=6)+ color_palette(palette = c("#577557", "#755757")) +
      theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust = 0.5))

We have an example in the DESeq2 workflow (rnaseqGene) where we have a heatmap of LFCs across time. It's the time course example.