Plotting change over time
1
0
Entering edit mode
mpuli011 • 0
@mpuli011-23216
Last seen 2.1 years ago

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)) deseq2 Plot • 348 views
0
Entering edit mode
@mikelove
Last seen 23 hours ago
United States

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

0
Entering edit mode

I see, so it would be a heat map as opposed to the other type of graph. Perfect, thank you. I have another question on a different topic, so I will ask it in another thread. Thank you.

0
Entering edit mode

I see, so it would be a heat map as opposed to the other type of graph. Perfect, thank you. I have another question on a different topic, so I will ask it in another thread. Thank you.

0
Entering edit mode

Hi again,

So I was able to create the heat map of change of taxa over time, but I noticed that the heat map is being created by using ASV ("OTU") ID. Unfortunately, it is not giving me the right information since I have multiple ID that are categorized under the same genus or order, etc. Is there any way for me to plot the heatmap retrieving the top 20 genus or taxa and then plotting that instead of the ASVID? I've looked through the forum and the outputs generated but I cannot seem to figure it out. Thanks

0
Entering edit mode

This seems like a question beyond the scope of DESeq2. I have to limit my time to just questions about usage of DESeq2.

0
Entering edit mode

I see, well can you tell me if I did the analysis correctly. I followed the tutorial you shared and when I change the name of the sequence to the actual name, there are replications of the taxa. When I graph by the sequence it makes sense because the ID are different, but if I wanted to plot by the gene name, per se, how can I do this.

Here is the graph w the edited information, where I replaced the sequence ID w the taxon name (in pdf).  0
Entering edit mode

No, i dont have time to review code.

0
Entering edit mode

Okay, well does Deseq cluster the gene ID's together if they are the same one and or not?

0
Entering edit mode

DESeq2 doesnt have any clustering functionality.

Do you mean pheatmap? This is a separate package. You should consult the help files for that package.